Taller 1 - Datos funcionales

Librerías

import os
import warnings
import numpy as np
import scipy.io as sio
import pandas as pd
import matplotlib.pyplot as plt
import skfda
from skfda.representation.basis import BSplineBasis
from skfda.preprocessing.smoothing import BasisSmoother
from skfda.exploratory.stats import mean, trim_mean, cov, depth_based_median
from skfda.exploratory.depth import ModifiedBandDepth, BandDepth
from tqdm import tqdm
from skfda.exploratory.visualization import Boxplot

# Suppress warnings
warnings.filterwarnings("ignore")

Cargar los datos

# Load MATLAB data
archivo_mat = "../Datos/data.mat"
datos = sio.loadmat(archivo_mat)

# Reshape data 
# Dimensions: 268 tests x 571 measurement points x 7 wavelengths
data = datos["X"].reshape(268, 571, 7, order='F')

# Wavelength labels
wavelengths = {
    0: 'Y1=230', 1: 'Y2=240', 2: 'Y3=255', 
    3: 'Y4=290', 4: 'Y5=305', 5: 'Y6=325', 6: 'Y7=340'
}

# Measurement points grid
dominio = np.arange(275, 560.5, 0.5)

# longitud de onda a trabajar
wavelength_idx = 6

Visualization of Raw Data

def plot_wavelength_data(organizado, n_samples=100):
    """
    Plot intensity data for first n_samples across 7 wavelengths
    """
    colors = plt.cm.tab10(np.linspace(0, 1, n_samples))
    fig, axes = plt.subplots(4, 2, figsize=(20, 25))
    fig.suptitle('7 Wavelengths - First 100 Samples', fontsize=16)
    
    axes = axes.ravel()
    x = np.arange(571)
    
    for wavelength in range(7):
        for sample in range(n_samples):
            data = organizado[sample, :, wavelength]
            axes[wavelength].plot(x, data, color=colors[sample], 
                                   linewidth=1.5, alpha=0.8)
        
        axes[wavelength].set_title(f'Wavelength {wavelength+1}', fontsize=12)
        axes[wavelength].set_xlabel('Measurement Points')
        axes[wavelength].set_ylabel('Intensity')
        axes[wavelength].grid(True, alpha=0.3)
    
    axes[-1].set_visible(False)
    plt.tight_layout(rect=[0, 0.03, 0.95, 0.95])
    plt.show()

plot_wavelength_data(data)

1.1 Suavizado spline

def smooth_functional_data(
    data, 
    wavelengths, 
    dominio, 
    base_path="../Resultados/punto_1", 
    n_basis=65, 
    best_param=1e-3, 
    rango_curvas=[None, None],
    mostrar_imagen_cada=None  
):
    """
    Perform functional data smoothing for specified range of samples and wavelengths
    
    Parameters:
    -----------
    data : numpy.ndarray
        Input data array with dimensions (n_samples, n_points, n_wavelengths)
    wavelengths : dict
        Dictionary of wavelengths
    dominio : array-like
        Domain points for the data
    base_path : str, optional
        Base directory to save results
    n_basis : int, optional
        Number of B-Spline basis functions
    best_param : float, optional
        Smoothing parameter (lambda)
    rango_curvas : list, optional
        Range of curves to process [start, end]. 
        If None, processes all curves. 
        Indices are 0-based.
    mostrar_imagen_cada : int or None, optional
        Show an image every 'n' processed samples.
        If None, no images are shown during processing.
    """
    # Crear la estructura de carpetas base
    if not os.path.exists(base_path):
        os.makedirs(base_path)

    # Crear subcarpetas para cada longitud de onda
    for wave_idx, wave_name in wavelengths.items():
        wave_path = os.path.join(base_path, wave_name)
        if not os.path.exists(wave_path):
            os.makedirs(wave_path)

    # Determinar el rango de curvas a procesar
    inicio = rango_curvas[0] if rango_curvas[0] is not None else 0
    fin = rango_curvas[1] if rango_curvas[1] is not None else 5  

    # Crear una barra de progreso
    total_iterations = fin - inicio
    with tqdm(total=total_iterations, desc="Procesando muestras") as pbar:
        # Contador de muestras procesadas
        muestras_procesadas = 0

        # Para cada longitud de onda
        for wave_idx, wave_name in wavelengths.items():
            # Si ya hemos procesado todas las muestras requeridas, salir del bucle
            if muestras_procesadas >= fin:
                break

            # Procesar muestras para esta longitud de onda
            for muestra in range(data.shape[0]):
                # Si ya hemos procesado todas las muestras requeridas, salir del bucle
                if muestras_procesadas >= fin:
                    break

                # Obtener los datos para esta muestra y longitud de onda
                X = data[muestra, :, wave_idx].flatten()
                
                # Crear una base B-Spline
                basis = BSplineBasis(domain_range=(275, 560), n_basis=n_basis)
                
                # Crear functional data grid
                fd = skfda.FDataGrid(data_matrix=X, grid_points=dominio)
                
                # Suavizado penalizado
                smoother = BasisSmoother(basis, smoothing_parameter=best_param)
                fd_smoothed = smoother.fit_transform(fd)
                
                # Decidir si mostrar la imagen
                if mostrar_imagen_cada is not None and muestras_procesadas % mostrar_imagen_cada == 0:
                    # Crear la figura
                    plt.figure(figsize=(12, 7))
                    
                    # Graficar los datos originales como puntos
                    plt.scatter(dominio, X, c='blue', s=30, alpha=0.5, label='Datos originales')
                    
                    # Graficar la curva ajustada
                    plt.plot(dominio, fd_smoothed(dominio)[0], 'r-',
                             label=f'Ajuste B-spline',
                             linewidth=2)
                    
                    # Personalización del gráfico
                    plt.title(f'Longitud de onda {wave_name} nm - Muestra global {muestra}\n'
                              f'Bases B-spline: {n_basis}', fontsize=12)
                    plt.xlabel('Longitud de onda (nm)', fontsize=11)
                    plt.ylabel('Intensidad', fontsize=11)
                    plt.legend(fontsize=10)
                    plt.grid(True, alpha=0.3)
                    
                    # Ajustar los márgenes
                    plt.tight_layout()
                    
                    # Mostrar la figura
                    plt.show()
                
                # Crear la figura para guardar
                plt.figure(figsize=(12, 7))
                
                # Graficar los datos originales como puntos
                plt.scatter(dominio, X, c='blue', s=30, alpha=0.5, label='Datos originales')
                
                # Graficar la curva ajustada
                plt.plot(dominio, fd_smoothed(dominio)[0], 'r-',
                         label=f'Ajuste B-spline',
                         linewidth=2)
                
                # Personalización del gráfico
                plt.title(f'Longitud de onda {wave_name} nm - Muestra global {muestra}\n'
                          f'Bases B-spline: {n_basis}', fontsize=12)
                plt.xlabel('Longitud de onda (nm)', fontsize=11)
                plt.ylabel('Intensidad', fontsize=11)
                plt.legend(fontsize=10)
                plt.grid(True, alpha=0.3)
                
                # Ajustar los márgenes
                plt.tight_layout()
                
                # Guardar la figura
                output_path = os.path.join(
                    base_path,
                    wave_name,
                    f'onda_{wave_idx}_muestra_global_{muestra}.png'
                )
                plt.savefig(output_path, dpi=300, bbox_inches='tight')
                plt.close()
                
                # Incrementar el contador de muestras procesadas
                muestras_procesadas += 1
                
                # Actualizar la barra de progreso
                pbar.update(1)

    print("Proceso completado. Se han generado todas las imágenes con el ajuste penalizado.")

total_curvas = 7*268
print(f'Hay un total de {total_curvas} curvas a procesar')
smooth_functional_data(data, wavelengths, dominio, rango_curvas=[0, 5], mostrar_imagen_cada=1)

Hay un total de 1876 curvas a procesar

Procesando muestras:   0%|          | 0/5 [00:00<?, ?it/s]

Procesando muestras:  20%|██        | 1/5 [00:01<00:06,  1.68s/it]

Procesando muestras:  40%|████      | 2/5 [00:03<00:05,  1.70s/it]

Procesando muestras:  60%|██████    | 3/5 [00:05<00:03,  1.72s/it]

Procesando muestras:  80%|████████  | 4/5 [00:07<00:01,  1.85s/it]

Procesando muestras: 100%|██████████| 5/5 [00:08<00:00,  1.77s/it]

Proceso completado. Se han generado todas las imágenes con el ajuste penalizado.

1.2

Estadistica descriptiva funcional

def plot_covariance_3d(cov_matrix):
    """
    Visualiza una matriz de covarianza en un gráfico 3D de superficie.
    
    Parámetros:
    cov_matrix (numpy.ndarray): Matriz de covarianza a visualizar
    """
    # Crear coordenadas para la malla
    x = np.arange(cov_matrix.shape[0])
    y = np.arange(cov_matrix.shape[1])
    X, Y = np.meshgrid(x, y)
    
    # Crear figura y eje 3D
    fig = plt.figure(figsize=(12, 8))
    ax = fig.add_subplot(111, projection='3d')
    
    # Graficar la superficie de la matriz de covarianza
    surf = ax.plot_surface(X, Y, cov_matrix, 
                           cmap='viridis',  # Puedes cambiar el colormap
                           edgecolor='none',
                           alpha=0.8)
    
    # Personalizar el gráfico
    ax.set_title('Visualización 3D de Matriz de Covarianza', fontsize=15)
    ax.set_xlabel('Dimensión 1', fontsize=12)
    ax.set_ylabel('Dimensión 2', fontsize=12)
    ax.set_zlabel('Valor de Covarianza', fontsize=12)
    
    # Añadir una barra de color
    fig.colorbar(surf, shrink=0.6, aspect=10)
    
    plt.tight_layout()
    plt.show()
    
def functional_stats_analysis(data, wavelength_idx):
    """
    Compute and visualize functional statistics
    """
    X = data[:, :, wavelength_idx]
    fd = skfda.FDataGrid(data_matrix=X, grid_points=dominio)
    
    # Mean computation.
    mean_func = mean(fd)
    mean_trim_func = trim_mean(fd, 0.1)
    
    # Plotting means
    plt.figure(figsize=(12, 6))
    for sample in range(268):
        plt.plot(dominio, X[sample], color="blue", linewidth=1.5, alpha=0.3)
    
    plt.plot(dominio, np.mean(X, axis=0), color='red', linewidth=3, label='Functional Mean')
    plt.plot(dominio, mean_trim_func(dominio)[0], color='green', linewidth=3, label='Trimmed Mean (10%)')
    
    plt.title(f'Individual Curves and Mean Functions - Wavelength {wavelengths[wavelength_idx]}')
    plt.xlabel('Wavelength (nm)')
    plt.ylabel('Intensity')
    plt.legend()
    plt.grid(True)
    plt.show()
    
    # Variance
    var_func = skfda.exploratory.stats.var(fd)
    var_func.dataset_name = f"Varianza Funcional"
    var_func.plot(color='purple', linewidth=2, legend=True)
    
    # Depth-based median
    median_func = depth_based_median(fd)
    median_func.dataset_name = f"Mediana basada en profundidad MBD"
    median_func.plot(color='orange', linewidth=2, legend=True)
    
    # covariance
    covariance = cov(fd, correction = 0)
    cov_matrix = covariance(fd.grid_points[0], fd.grid_points[0])
    plot_covariance_3d(cov_matrix)
    return fd

fd_analysis = functional_stats_analysis(data, wavelength_idx)

Profundidad y class_outliners

Analisis de profundidad funcional por BD y MBD

Definición de Band Depth (BD)

La Band Depth (BD) mide la "centralidad" de una función $x$ en relación con un conjunto de curvas $x_1, x_2, \ldots, x_n$. Se basa en calcular la proporción de bandas, formadas por subconjuntos de $j$ curvas, que contienen completamente a la función $x$.

1. Bandas Formadas por Curvas: Se tiene que el gráfico de la función $x$ es el subconjunto del plano definido por $G(x) = \{(t, x(t)) : t \in I\}.$ Para $j \geq 2$, una banda $B(x_{i_1}, x_{i_2}, \ldots, x_{i_j})$ delimitada por $j$ curvas $x_{i_1}, \ldots, x_{i_j}$ se define como:
   $$ B(x_{i_1}, x_{i_2}, \ldots, x_{i_j}) = \{ (t, y) : t \in I, \min_{r=1, \ldots, j} x_{i_r}(t) \leq y \leq \max_{r=1, \ldots, j} x_{i_r}(t) \}. $$

2. Proporción de Inclusión:
   El índice de profundidad basado en bandas para $x$ es:
   $$ S_n^{(j)}(x) = \binom{n}{j}^{-1} \sum_{1 \leq i_1 < i_2 < \cdots < i_j \leq n} I\{ G(x) \subseteq B(x_{i_1}, x_{i_2}, \ldots, x_{i_j}) \}, $$
   donde $I\{\cdot\}$ es la función indicadora que vale 1 si $G(x)$ (la gráfica de $x$) está contenida en la banda, y 0 en caso contrario. Observe que la cantidad anterior calcula la proporción de bandas que contienen a la curva $x$.

3. Profundidad:
   La profundidad de banda para $J \geq 2$ se define como:
   $$ S_{n,J}(x) = \sum_{j=2}^{J} S_n^{(j)}(x). $$

Modified Band Depth (MBD)

La Modified Band Depth (MBD) extiende la definición de BD al permitir medir qué proporción del intervalo $I$ está contenida dentro de una banda. Esto hace que la métrica sea más robusta a comportamientos locales irregulares.

1. Conjunto de Inclusión Parcial:
   Se define el conjunto donde la función $x$ está dentro de la banda como:
   $$ A_j(x; x_{i_1}, \ldots, x_{i_j}) = \{ t \in I : \min_{r=1, \ldots, j} x_{i_r}(t) \leq x(t) \leq \max_{r=1, \ldots, j} x_{i_r}(t) \}. $$

2. Medida Proporcional:
   Se mide la proporción de tiempo que $x$ está dentro de la banda:
   $$ \lambda_r(A_j(x)) = \frac{\lambda(A_j(x))}{\lambda(I)}, $$
   donde $\lambda$ es la medida de Lebesgue.

3. Profundidad Generalizada:
   La MBD es:
   $$ GS_n^{(j)}(x) = \binom{n}{j}^{-1} \sum_{1 \leq i_1 < \cdots < i_j \leq n} \lambda_r(A_j(x; x_{i_1}, \ldots, x_{i_j})). $$
   Finalmente:
   $$ GS_{n,J}(x) = \sum_{j=2}^{J} GS_n^{(j)}(x). $$

Implementación del Functional Boxplot

1. Ordenación por Profundidad:
   Ordena las curvas de acuerdo con su profundidad BD o MBD, desde la más profunda hasta la menos profunda.

2. Mediana Funcional:
   La curva más profunda se define como la mediana funcional:
   $$ \hat{m}_n = \arg\max_{x \in \{x_1, \ldots, x_n\}} S_n^{(J)}(x) \quad \text{(o MBD)}. $$

3. Región Central y Outliers:

   • Define una banda central basada en el rango intercuartílico de las curvas más profundas.
   • Las curvas fuera de esta banda se consideran outliers.

4. Visualización:

   • La mediana funcional y las bandas centrales se visualizan de manera análoga al boxplot estándar, representando las curvas más centrales y las más extremas.

Propiedades y Ventajas

• Robustez: Ambas métricas son resistentes a valores atípicos, pero la MBD es más flexible frente a irregularidades locales.

• Eficiencia Computacional: La BD y MBD son rápidas de calcular, incluso para datos de alta dimensión.

• Aplicabilidad: Son útiles para detectar outliers y analizar la centralidad de curvas en conjuntos funcionales.

• Implementación: Se encuentra implementada en librerias de R o Python, permitiendo calcular estas métricas y visualizar el functional boxplot.

def agregar_valores(data, values, wavelength):
    """
    Agrega un nuevo valor al final de cada fila en la dimensión de wavelength.
    
    Parámetros:
    - data: Array numpy con dimensiones (268, 571, ...)
    - values: Array numpy con dimensiones (268,)
    - wavelength: Índice de la longitud de onda donde se agregará el nuevo valor
    
    Retorna:
    - Nuevo array numpy con valores agregados
    """
    # Crear un nuevo array con un espacio adicional en la dimensión de wavelength
    data_augm = np.zeros((data.shape[0], 
                                      data.shape[1] + 1, 
                                      data.shape[2]))
    
    # Verificar dimensiones
    if values.shape[0] != data.shape[0]:
        raise ValueError("Las dimensiones de 'values' y 'data' no coinciden")
    
    # Copiar los datos originales
    data_augm[:, :data.shape[1], :] = data
    
    # Agregar los valores al final
    for i in range(data.shape[0]):
        data_augm[i, -1, wavelength] = values[i]
    
    return data_augm

def analisis_profundidad_funcional(data, fd, wavelength, metodo='mbd'):
    """
    Realiza análisis de profundidad funcional y genera visualización de percentiles.
    
    Parámetros:
    - data: Array numpy con dimensiones (268, 571, ...)
    - fd: Datos funcionales para cálculo de profundidad
    - wavelength: Índice de la longitud de onda 
    - metodo: 'mbd' (ModifiedBandDepth) o 'bd' (BandDepth)
    
    Retorna:
    - DataFrame con valores de profundidad ordenados
    - Gráfico de percentiles
    """
    # Selección del método de profundidad
    if metodo == 'mbd':
        depth = skfda.exploratory.depth.ModifiedBandDepth()
    elif metodo == 'bd':
        depth = skfda.exploratory.depth.BandDepth()
    else:
        raise ValueError("Método de profundidad no válido. Use 'mbd' o 'bd'")
    
    # Cálculo de valores de profundidad
    values = depth(fd)
    
    # Agregar valores de profundidad al array original
    data_modificado = agregar_valores(data, values, wavelength)
    
    # Crear DataFrame
    df = pd.DataFrame(data_modificado[:,:, wavelength])
    df.sort_values(by=571, axis=0, ascending=False, inplace=True)
    
    # imprime las 10 primeras filas de df
    print(f"Valores de profundidad ordenados - Método {metodo.capitalize()}")
    display(df.head(20))
    
    
    # Cálculo de percentiles
    percentile_90 = df.quantile(0.90)[:570]
    percentile_95 = df.quantile(0.95)[:570]
    
    # Gráfico de percentiles
    plt.figure(figsize=(12, 6))
    plt.plot(percentile_90, label='Percentil 90%', color='blue', linestyle='--')
    plt.plot(percentile_95, label='Percentil 95%', color='red', linestyle='--')
    plt.title(f'Percentiles 90% y 95% - Profundidad {metodo.capitalize()}')
    plt.xlabel('Muestras')
    plt.ylabel('Profundidad')
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()
    
    return df

df_bd = analisis_profundidad_funcional(data, fd_analysis, wavelength_idx, metodo='bd')
df_mbd = analisis_profundidad_funcional(data, fd_analysis, wavelength_idx, metodo='mbd')

Valores de profundidad ordenados - Método Bd

	0	1	2	3	4	5	6	7	8	9	...	562	563	564	565	566	567	568	569	570	571
26	1.906	2.055	2.294	2.658	3.165	3.816	4.624	5.627	6.872	8.397	...	16.591	16.396	16.286	16.302	16.418	16.553	16.668	16.818	17.069	0.046816
57	1.627	1.808	2.112	2.540	3.078	3.736	4.568	5.617	6.853	8.180	...	16.564	16.182	15.986	16.019	16.187	16.410	16.766	17.419	18.385	0.042680
64	1.891	1.959	2.183	2.578	3.106	3.735	4.483	5.409	6.565	7.970	...	16.935	16.715	16.780	17.000	17.185	17.216	17.104	16.953	16.866	0.041897
123	1.660	1.886	2.191	2.551	2.971	3.490	4.156	4.997	6.023	7.219	...	17.802	17.634	17.488	17.346	17.277	17.372	17.643	18.008	18.363	0.039298
31	1.932	2.032	2.262	2.628	3.111	3.716	4.484	5.459	6.659	8.093	...	16.733	16.472	16.290	16.276	16.447	16.756	17.160	17.656	18.218	0.038376
75	1.995	2.177	2.438	2.816	3.355	4.073	4.958	5.991	7.184	8.602	...	16.807	16.785	16.790	16.789	16.796	16.894	17.144	17.467	17.694	0.037677
36	1.940	2.146	2.480	2.939	3.506	4.160	4.918	5.853	7.059	8.554	...	17.275	17.381	17.550	17.657	17.623	17.511	17.501	17.740	18.218	0.036140
93	1.697	1.802	2.037	2.414	2.909	3.517	4.286	5.269	6.474	7.889	...	17.627	17.714	17.588	17.368	17.186	17.111	17.178	17.437	17.893	0.035217
74	1.872	1.965	2.206	2.622	3.193	3.887	4.725	5.796	7.200	8.970	...	16.162	16.119	16.264	16.506	16.724	16.847	16.880	16.900	16.983	0.034714
59	1.746	1.968	2.299	2.721	3.209	3.767	4.440	5.291	6.367	7.691	...	16.393	16.341	16.243	16.101	15.934	15.793	15.775	15.966	16.352	0.034546
79	1.859	2.115	2.443	2.845	3.349	3.997	4.853	5.999	7.487	9.277	...	17.132	17.309	17.461	17.541	17.591	17.698	17.946	18.392	19.004	0.033876
98	1.642	1.846	2.138	2.532	3.057	3.745	4.602	5.635	6.871	8.334	...	16.378	16.411	16.468	16.487	16.494	16.537	16.616	16.702	16.775	0.033876
32	1.648	1.797	2.060	2.470	3.030	3.702	4.453	5.299	6.303	7.528	...	16.816	16.854	16.888	16.875	16.862	16.961	17.251	17.731	18.323	0.033177
126	1.654	1.893	2.263	2.738	3.276	3.854	4.492	5.245	6.181	7.357	...	17.201	16.850	16.612	16.590	16.720	16.870	16.992	17.149	17.406	0.033121
52	1.743	1.888	2.183	2.645	3.234	3.889	4.607	5.478	6.627	8.100	...	16.137	16.129	16.054	15.941	15.900	15.985	16.047	15.816	15.179	0.032730
35	2.002	2.110	2.362	2.819	3.474	4.247	5.078	5.987	7.047	8.330	...	16.389	16.181	16.098	16.166	16.403	16.771	17.128	17.328	17.353	0.032562
41	1.978	2.146	2.425	2.849	3.412	4.094	4.904	5.885	7.073	8.482	...	17.442	17.466	17.420	17.310	17.169	17.051	16.997	17.009	17.053	0.032394
61	1.933	2.077	2.295	2.618	3.108	3.831	4.806	5.980	7.266	8.621	...	17.420	17.296	17.296	17.389	17.489	17.588	17.779	18.159	18.717	0.032282
34	1.958	2.099	2.325	2.683	3.209	3.897	4.719	5.690	6.880	8.334	...	17.422	17.298	17.194	17.179	17.250	17.391	17.624	17.962	18.362	0.032282
55	1.737	1.982	2.370	2.903	3.561	4.327	5.205	6.205	7.339	8.652	...	17.794	17.599	17.474	17.416	17.453	17.621	17.890	18.156	18.334	0.032059

20 rows × 572 columns

Valores de profundidad ordenados - Método Mbd

	0	1	2	3	4	5	6	7	8	9	...	562	563	564	565	566	567	568	569	570	571
174	1.362	1.439	1.614	1.876	2.189	2.551	3.033	3.726	4.668	5.815	...	15.712	15.364	15.146	15.063	15.070	15.139	15.312	15.672	16.221	0.484090
144	1.462	1.797	1.863	2.483	2.945	3.007	3.613	4.068	4.912	6.689	...	15.652	16.044	16.044	16.074	17.268	17.211	15.376	15.672	17.251	0.478247
78	2.024	2.117	2.277	2.538	2.941	3.516	4.269	5.212	6.386	7.833	...	16.231	16.252	16.280	16.320	16.392	16.519	16.751	17.144	17.677	0.477717
183	1.365	1.445	1.598	1.824	2.124	2.527	3.080	3.813	4.704	5.695	...	15.607	15.488	15.427	15.448	15.538	15.661	15.772	15.848	15.905	0.477697
96	2.079	2.198	2.412	2.761	3.262	3.892	4.627	5.475	6.482	7.715	...	16.176	16.130	16.258	16.501	16.750	16.876	16.773	16.424	15.916	0.476309
43	1.734	1.904	2.129	2.426	2.825	3.367	4.093	5.035	6.202	7.575	...	16.413	16.335	16.276	16.220	16.188	16.240	16.438	16.804	17.301	0.475230
52	1.743	1.888	2.183	2.645	3.234	3.889	4.607	5.478	6.627	8.100	...	16.137	16.129	16.054	15.941	15.900	15.985	16.047	15.816	15.179	0.475000
138	1.529	1.667	1.918	2.315	2.861	3.527	4.268	5.067	5.944	6.945	...	16.178	16.122	16.134	16.181	16.187	16.185	16.331	16.759	17.422	0.472879
27	1.824	2.006	2.313	2.709	3.163	3.705	4.405	5.328	6.501	7.911	...	16.060	16.091	16.072	16.111	16.263	16.425	16.424	16.198	15.850	0.470869
175	1.399	1.476	1.632	1.859	2.154	2.533	3.029	3.665	4.441	5.342	...	16.526	16.405	16.435	16.556	16.617	16.541	16.430	16.450	16.646	0.467707
180	1.989	1.882	1.906	2.129	2.540	3.075	3.680	4.358	5.155	6.111	...	16.591	16.528	16.421	16.276	16.125	16.002	15.948	15.993	16.126	0.467657
65	2.115	2.065	2.121	2.358	2.797	3.392	4.093	4.912	5.915	7.143	...	16.268	16.199	16.073	15.990	16.019	16.114	16.116	15.866	15.351	0.467653
25	2.122	2.153	2.351	2.730	3.247	3.876	4.661	5.678	6.952	8.438	...	16.246	15.890	15.689	15.721	15.980	16.356	16.674	16.796	16.705	0.467151
158	1.452	1.411	1.464	1.722	2.169	3.570	2.911	3.319	5.210	6.077	...	16.954	16.984	17.875	15.328	16.087	16.912	15.874	16.366	17.390	0.466007
22	1.836	2.036	2.395	2.882	3.454	4.125	4.959	5.988	7.185	8.527	...	17.130	17.183	17.113	17.008	16.953	16.958	17.018	17.161	17.398	0.464580
100	2.249	2.330	2.561	2.930	3.398	3.954	4.634	5.509	6.656	8.121	...	16.295	16.342	16.505	16.708	16.874	16.970	17.041	17.168	17.381	0.463610
98	1.642	1.846	2.138	2.532	3.057	3.745	4.602	5.635	6.871	8.334	...	16.378	16.411	16.468	16.487	16.494	16.537	16.616	16.702	16.775	0.463585
181	1.269	1.508	1.783	2.075	2.399	2.791	3.273	3.853	4.547	5.406	...	16.897	16.397	16.099	16.144	16.434	16.714	16.705	16.237	15.372	0.462293
143	1.514	1.655	1.883	2.220	2.676	3.259	3.979	4.830	5.777	6.799	...	16.120	15.983	15.836	15.721	15.659	15.616	15.503	15.260	14.920	0.461448
178	1.723	1.638	1.620	1.738	2.020	2.439	2.950	3.554	4.295	5.205	...	16.595	16.499	16.396	16.349	16.355	16.338	16.265	16.219	16.294	0.460649

20 rows × 572 columns

Boxplot funcional

Construcción de Functional Boxplots

Un functional boxplot extiende el concepto del boxplot clásico a datos funcionales y consta de los siguientes elementos:

1. Región central del 50%: Representa el rango central de las curvas más profundas, delimitada por la banda que contiene el 50% de las curvas más centrales.
   $$ C_{0.5} = \{ (t, y(t)) : \min_{r=1,\ldots,\lfloor n/2 \rfloor} y_{[r]}(t) \leq y(t) \leq \max_{r=1,\ldots,\lfloor n/2 \rfloor} y_{[r]}(t) \}. $$

2. Curva mediana funcional: La curva más central basada en la profundidad (por ejemplo, BD o MBD).
   $$ \hat{m}_n = \arg\max_{y \in \{y_1, \ldots, y_n\}} \text{BD}_n(y). $$

3. Región máxima no atípica: El rango completo de curvas no consideradas outliers.

4. Valores atípicos: Curvas fuera de los límites establecidos por la región central ampliada en 1.5 veces su rango.

Procedimiento de Construcción

1. Ordenación: Las curvas se ordenan según su profundidad (BD o MBD) de forma descendente.

2. Cálculo de la región central: Se identifica la banda que contiene las curvas más profundas (usualmente el 50%).

3. Detección de outliers: Se identifican curvas fuera de los límites ampliados.

4. Visualización:

   • La región central del 50% se representa como el análogo funcional del rango intercuartílico (IQR).
   • La mediana funcional se muestra como la curva más representativa.
   • Los valores atípicos se destacan.

Conclusiones

• Robustez: El functional boxplot es robusto frente a valores atípicos y proporciona una representación visual confiable de la dispersión y centralidad de las curvas.
• Descriptivo y Comparativo: Permite comparar poblaciones de curvas y detectar diferencias significativas entre conjuntos de datos funcionales.
• Aplicaciones Versátiles: Útil en múltiples campos como datos temporales, espaciales y espaciotemporales.

def functional_boxplot(data, wavelength_idx, F=1.5):
    """
    Create functional boxplot with outliers
    """
    
    X = data[:, :, wavelength_idx]
    fd = skfda.FDataGrid(data_matrix=X, grid_points=dominio, dataset_name=f'Boxplot with F: {F}')
    fdBoxplot = Boxplot(fd, factor=F)        
    fdBoxplot.plot()
    return fdBoxplot
    

fboxplot = functional_boxplot(data, wavelength_idx)
outliers = fboxplot.outliers
# Compara los outliers con data y muestrame los valores de data que son outliers
data_outliers = pd.DataFrame(data[outliers, :, wavelength_idx])
display(data_outliers)

	0	1	2	3	4	5	6	7	8	9	...	561	562	563	564	565	566	567	568	569	570
0	1.891	2.045	2.286	2.603	2.995	3.497	4.161	5.039	6.181	7.640	...	25.238	25.418	25.498	25.481	25.450	25.457	25.492	25.541	25.645	25.852
1	1.966	2.171	2.497	2.931	3.456	4.090	4.863	5.796	6.910	8.263	...	22.321	22.459	22.633	22.694	22.543	22.194	21.763	21.392	21.164	21.059
2	1.206	1.321	1.466	1.650	1.884	2.183	2.580	3.117	3.816	4.655	...	24.670	24.268	24.006	23.916	23.959	24.051	24.130	24.230	24.465	24.887

3 rows × 571 columns

Boxplot funcional adjustado

Ajuste del Factor en Functional Boxplots

El factor constante 1.5, utilizado tradicionalmente en boxplots clásicos para la detección de outliers, también se aplica en functional boxplots como multiplicador del rango de la región central del 50%. Sin embargo, este valor fijo puede no ser adecuado en todos los contextos debido a la naturaleza de los datos funcionales.

Razones para Ajustar el Factor

1. Dependencia Temporal y Espaciotemporal:

   • Las curvas funcionales suelen presentar dependencia temporal (en datos univariados) y espaciotemporal (en datos multivariados).
   • Esta dependencia puede alterar la dispersión observada en las curvas y, por tanto, afecta la detección de outliers si se usa un factor fijo como 1.5.

2. Distribución de los Outliers:

   • En un boxplot clásico, el factor 1.5 está basado en la distribución normal estándar, donde la probabilidad de detectar un valor como outlier es aproximadamente 0.7%.
   • Este factor no toma en cuenta la correlación inherente en los datos funcionales, lo que puede llevar a sobreestimación o subestimación de outliers.

3. Impacto en la Visualización y Detección:

   • Un factor fijo puede no ser adecuado si las curvas tienen una estructura de correlación compleja.
   • Ajustar el factor permite que la probabilidad de detectar valores atípicos se adapte mejor a las características del conjunto de datos.

Propuesta de Ajuste

El factor puede ajustarse dinámicamente para garantizar que la probabilidad de no detectar outliers sea consistente con la estructura de los datos. Por ejemplo:

• En ausencia de outliers reales, ajustar el factor para que la probabilidad de no detectar outliers sea del 99.3%.
• Este enfoque mejora la sensibilidad del functional boxplot para identificar valores atípicos mientras mantiene la robustez frente a la dependencia.

Simulación

Simularemos curvas, a apartir de un proceso Gaussiano, definiendo la estructura de variación a partir de una matriz de correlaciones basadas en un estimador robusto de la varianza y covarianza dado por la siguiente espresión:

Estimador robusto de una matriz de dispersión (Ma & Genton, 2001)

El estimador de covarianza robusta está dado por:

$$ \text{Cov}(X, Y) = \frac{\sigma_X \sigma_Y}{4} \left[ \text{Var}\left( \frac{X}{\sigma_X} + \frac{Y}{\sigma_Y} \right) - \text{Var}\left( \frac{X}{\sigma_X} - \frac{Y}{\sigma_Y} \right) \right], $$

Este estimador reduce la influencia de outliers en la matriz de dispersión, crucial en escenarios con dependencia estructural en los datos.

Estimador altamente robusto de escala (Rousseeuw y Croux, 1992, 1993)

Dado un conjunto de datos $\{z_1, \ldots, z_n\}$, el estimador robusto de escala es:

$$ \hat{\sigma}_Z = Q_n(Z) = d \, \Big[ \lvert z_i - z_j \rvert; \, i < j, \, i, j = 1, 2, \ldots, n \Big]_{(k)}, $$

donde:

• $d = 2.2191$ es un factor de consistencia,
• $k = \lfloor ((n/2) + 2)/4 \rfloor + 1$,
• $(k)$ indica aproximadamente el primer cuartil para muestras grandes.

Este método permite medir la escala de los datos mientras mantiene robustez frente a valores atípicos.

Una vez que se han calculado la matriz $Cov(X,Y)$ con la data inicial, se procede a hacer una serie de simulaciones para encontrar el mejor F. Para esto simulamos 1000 matrices de muestras de los datos, y calculamos

Conclusión

El ajuste del factor en functional boxplots es crucial para asegurar una detección de outliers confiable en escenarios con alta dependencia temporal o espaciotemporal. Aunque el valor 1.5 funciona bien en contextos simples, en datos funcionales complejos, un factor adaptable proporciona mejores resultados sin comprometer la visualización.

import numpy as np
import skfda
from skfda.exploratory.outliers import BoxplotOutlierDetector
import numpy as np
import scipy.special as sc


def calculate_stats(data):
    """
    Calcula un cuantil específico de las diferencias absolutas para cada columna
    de una matriz 2D, considerando únicamente los pares donde i < j.

    Args:
        data (numpy.ndarray): Una matriz 2D donde cada columna representa una variable.

    Returns:
        numpy.ndarray: Un arreglo con los cuantiles calculados para cada columna.

    Raises:
        ValueError: Si la entrada no es una matriz 2D numérica.
    """
    if not isinstance(data, np.ndarray) or data.ndim != 2:
        raise ValueError("Input data must be a 2D NumPy array.")
    if not np.issubdtype(data.dtype, np.number):
        raise ValueError("Input data must be numeric.")

    n_cols = data.shape[1]
    d = 2.2191
    quantiles = []
    

    for col_idx in range(n_cols):
        col_data = data[:, col_idx]
        n = len(col_data)
        diffs = []

        # Calcular diferencias solo para i < j
        for i in range(n):
            for j in range(i + 1, n):  # Solo i < j
                diffs.append(np.abs(col_data[i] - col_data[j]))

        diffs = np.array(diffs)
        # Calcular el índice del cuantil específico
        k = int(np.floor((sc.comb(n, 2) + 2) / 4)) + 1
        # multiplica k por d
        k_int = int(k * d)
        quantiles.append(np.partition(diffs, k - 1)[k - 1])  # Obtener el k-ésimo valor más pequeño

    return np.array(quantiles)


def calculate_covariance_matrix(array_reshaped, stats):
    # Número de columnas
    num_columns = array_reshaped.shape[1]

    # Inicializar la matriz de covarianza
    covariance_matrix = np.zeros((num_columns, num_columns))

    for i in range(num_columns):
        for j in range(num_columns):
            # Extraer las columnas i y j
            X = array_reshaped[:, i, 4]
            Y = array_reshaped[:, j, 4]

            # Obtener sigma_X y sigma_Y de stats
            sigma_X = stats[i]
            sigma_Y = stats[j]

            # Calcular términos de la fórmula
            term_1 = (X / sigma_X) + (Y / sigma_Y)
            term_2 = (X / sigma_X) - (Y / sigma_Y)

            # Calcular las varianzas
            var_term_1 = np.var(term_1)
            var_term_2 = np.var(term_2)

            # Aplicar la fórmula
            covariance_matrix[i, j] = (sigma_X * sigma_Y / 4) * (var_term_1 - var_term_2)

    return covariance_matrix




def calculate_outlier_proportions(covariance_matrix, num_matrices=1000, num_samples=100, num_groups=10):
    # Generar las 1000 matrices de muestras
    concatenated_samples = [
        np.random.multivariate_normal(np.zeros(covariance_matrix.shape[0]), covariance_matrix, size=num_samples)
        for _ in range(num_matrices)
    ]

    # Dividir las matrices en 10 grupos
    matrices_per_group = num_matrices // num_groups
    groups = [
        concatenated_samples[i * matrices_per_group:(i + 1) * matrices_per_group]
        for i in range(num_groups)
    ]

    # Configurar los valores de factor para cada grupo
    factors = [0.9 + 0.1 * i for i in range(num_groups)]

    # Inicializar resultados
    group_proportions = []

    # Procesar cada grupo de matrices
    for group_idx, (group, factor) in enumerate(zip(groups, factors)):
        proportions = []

        for sample in group:
            # Crear el objeto FDataGrid para la matriz actual
            num_columns = covariance_matrix.shape[0]
            grid_points = np.arange(275, 275 + 0.5 * num_columns, 0.5)
            fd = skfda.FDataGrid(data_matrix=sample, grid_points=grid_points)

            # Aplicar el detector de outliers
            out_detector = BoxplotOutlierDetector(factor=factor)
            outlier_predictions = out_detector.fit_predict(fd)

            # Calcular proporción de outliers
            proportion = np.sum(outlier_predictions == -1) / len(outlier_predictions)
            proportions.append(proportion)

        # Calcular promedio de proporciones para el grupo actual
        group_avg_proportion = 1-np.mean(proportions)
        group_proportions.append(group_avg_proportion)
        print(f"Grupo {group_idx + 1} con factor {factor:.1f}: Promedio de proporciones de outliers = {group_avg_proportion:.4f}")

    return group_proportions

stats = calculate_stats(data[:,:,wavelength_idx])
covariance_matrix = calculate_covariance_matrix(data, stats)
calculate_outlier_proportions(covariance_matrix, num_matrices=1000, num_samples=100, num_groups=10)

Grupo 1 con factor 0.9: Promedio de proporciones de outliers = 0.9603
Grupo 2 con factor 1.0: Promedio de proporciones de outliers = 0.9659
Grupo 3 con factor 1.1: Promedio de proporciones de outliers = 0.9782
Grupo 4 con factor 1.2: Promedio de proporciones de outliers = 0.9878
Grupo 5 con factor 1.3: Promedio de proporciones de outliers = 0.9906
Grupo 6 con factor 1.4: Promedio de proporciones de outliers = 0.9941
Grupo 7 con factor 1.5: Promedio de proporciones de outliers = 0.9976
Grupo 8 con factor 1.6: Promedio de proporciones de outliers = 0.9983
Grupo 9 con factor 1.7: Promedio de proporciones de outliers = 0.9989
Grupo 10 con factor 1.8: Promedio de proporciones de outliers = 0.9989

[np.float64(0.9603),
 np.float64(0.9659),
 np.float64(0.9782),
 np.float64(0.9878),
 np.float64(0.9906),
 np.float64(0.9941),
 np.float64(0.9976),
 np.float64(0.9983),
 np.float64(0.9989),
 np.float64(0.9989)]

fboxplot = functional_boxplot(data, wavelength_idx, F=1.4)
outliers = fboxplot.outliers
# Compara los outliers con data y muestrame los valores de data que son outliers
data_outliers = pd.DataFrame(data[outliers, :, wavelength_idx])

Análisis de la Profundidad de Variación Total y Similitud de Forma

** 1. Introducción**

Sea $X$ un proceso estocástico en $\tau$, donde $\tau$ es un intervalo en $\mathbb{R}$, con distribución $F_X$. Denotamos:

• $f$: una función.
• $f(t)$: el valor funcional en un tiempo $t$ dado.

2. Profundidad de Variación Total Puntual

Para una función $f(t)$, definimos:

$$ R_f(t) = \mathbb{I}_{\{ X(t) \leq f(t) \}}, $$

donde $\mathbb{I}$ es la función indicadora. A partir de esta relación:

$$ p_f(t) = \mathbb{E}[R_f(t)] = \mathbb{P}\{ X(t) \leq f(t) \}, $$

lo que indica la posición de $f(t)$ con respecto a $X(t)$.

La profundidad de variación total puntual se define como:

$$ D_f(t) = \text{Var}[R_f(t)] = p_f(t)(1 - p_f(t)). $$

3. Descomposición de la Profundidad de Variación Total

Sea $s, t$ dos puntos en el tiempo donde $s = t - \Delta$. La profundidad de variación total puntual tiene la siguiente descomposición:

$$ D_f(t) = \text{Var}[R_f(t)] = \text{Var}[\mathbb{E}\{R_f(t) \mid R_f(s)\}] + \mathbb{E}[\text{Var}\{R_f(t) \mid R_f(s)\}]. $$

Esta descomposición implica que la varianza total de $R_f(t)$ se puede dividir en:

• $\text{Var}[\mathbb{E}\{R_f(t) \mid R_f(s)\}]$: la componente de forma, que representa la porción de variabilidad de $R_f(t)$ explicada por $R_f(s)$.
• $\mathbb{E}[\text{Var}\{R_f(t) \mid R_f(s)\}]$: la componente de magnitud, que es independiente de $R_f(s)$.

4. Profundidad de Variación Total Funcional (TVD)

La profundidad de variación total funcional (TVD, por sus siglas en inglés) para una función $f$ es:

$$ TVD(f) = \int w(t) D_f(t) \, dt, $$

donde $w(t)$ es una función de peso definida sobre el dominio de $f(t)$. Si consideramos que $w(t)$ es constante, tal que:

$$ w(t) \equiv \frac{1}{|\tau|}, $$

entonces el TVD se comporta de manera similar a la profundidad de banda modificada.

5. Similitud de Forma (Shape Similarity, SS)

Para cualquier función $f$ en $\mathbb{R}$, la similitud de forma dada una ventana temporal $\Delta$ se define como:

$$ SS(f; \Delta) = \int v(t; \Delta) S_f(t; \Delta) \, dt, $$

donde:

• $v(t; \Delta)$: función de peso, dada por:

$$ v(t; \Delta) = \frac{|f(t) - f(t - \Delta)|}{\int |f(t) - f(t - \Delta)|}. $$

• $S_f(t; \Delta)$: componente de forma, definida como:

$$ S_f(t; \Delta) = \begin{cases} \frac{\text{Var}[R_f(t) \mid R_f(t - \Delta)]}{D_f(t)}, & \text{si } D_f(t) \neq 0, \\ 1, & \text{si } D_f(t) = 0. \end{cases} $$

Interpretación de la Similitud de Forma

Valores pequeños de $SS(f; \Delta)$ están asociados con una mayor desviación en la forma. Sin embargo, para pares atípicos $(f(t - \Delta), f(t))$, si el denominador $D_f(t)$ es muy pequeño, $S_f(t; \Delta)$ podría no reflejar adecuadamente la desviación.

Para capturar mejor esta desviación, se centra el par $(f(t - \Delta), f(t))$ mediante un desplazamiento $\delta_t$, donde:

$$ \delta_t = f(t) - \text{mediana}\{X(t)\}. $$

El par ajustado es:

$$ (f(t - \Delta), f(t)) - \delta_t. $$

6. Similitud de Forma Modificada (Modified Shape Similarity, MSS)

La similitud de forma modificada (MSS) para una función $f$, dada una ventana temporal fija $\Delta$, se define como:

$$ MSS(f; \Delta) = \int v(t; \Delta) S_f(t; \Delta) \, dt, $$

donde $S_f(t; \Delta)$ es:

$$ S_f(t; \Delta) = \frac{\text{Var}[R_f(t) \mid R_f(t - \Delta)]}{D_f(t)}. $$

Además, se define $f(s; \Delta)$ como:

$$ f(s; \Delta) = \begin{cases} f(s) - f(s + \Delta) + \text{mediana}\{X(s + \Delta)\}, & \text{si } s \in \mathbb{R}. \end{cases} $$

---

Outlier Detection Procedure

The process for detecting outliers among a set of $n$ sample curves involves the following steps:

1. Compute Metrics

   • Estimate the total variation depth and modified shape similarity for each curve as before

2. Identify Shape Outliers

   • Create a classical boxplot for the $n$ values of the modified shape similarity.
   • Detect outliers using the $F \times IQR$ empirical rule, where $F$ is a user-adjustable factor.
   • Curves with modified shape similarity values below the lower fence are classified as shape outliers.

3. Identify Magnitude Outliers

   • Remove shape outliers and construct a functional boxplot using the total variation depth.
   • Detect magnitude outliers by identifying curves outside the 50% central region (based on the original observations) expanded by a factor of 1.5. This factor can be adjusted using bootstrap methods .

import rpy2

%load_ext rpy2.ipython

X = data[:, :, wavelength_idx]

%R -i X

Además: Aviso:
In (function (package, help, pos = 2, lib.loc = NULL, character.only = FALSE,  :
  library ‘/usr/lib/R/site-library’ contains no packages

%%R
library("fdaoutlier")
# Just show me the first 5 rows in R
class_outliners <- tvdmss(X)

%R -o class_outliners

def plot_kind_class_outliners(outliers, X):
    """
    Plot outliers and non-outliers with different colors and labels
    """
    # Convierte los índices de outliers a un conjunto para una búsqueda eficiente
    magnitude_outliers_set = set(outliers['magnitude_outliers'])
    shape_outliers_set = set(outliers['shape_outliers'])
    
    # Crear máscaras para diferentes tipos de outliers
    magnitude_outliers_mask = np.array([i in magnitude_outliers_set for i in range(len(X))])
    shape_outliers_mask = np.array([i in shape_outliers_set for i in range(len(X))])
    
    # Encuentra los magnitude outliers en X
    magnitude_outliers = X[magnitude_outliers_mask, :]
    
    # Encuentra los shape outliers en X
    shape_outliers = X[shape_outliers_mask, :]
    
    # Encuentra las curvas que no son outliers
    non_outliers_mask = ~(magnitude_outliers_mask | shape_outliers_mask)
    non_outliers = X[non_outliers_mask, :]
    
    # Grafica los tipos de curvas
    plt.figure(figsize=(12, 6))
    
        
    for i in range(shape_outliers.shape[0]):
        plt.plot(shape_outliers[i], color='blue', alpha=0.7, label=('' if i==0 else '_') + 'Outliers de Forma')
    
    for i in range(magnitude_outliers.shape[0]):
        plt.plot(magnitude_outliers[i], color='red', alpha=0.7, label=('' if i==0 else '_') + 'Outliers de Magnitud')
        
    plt.title('Outliers de Magnitud y Forma')
    plt.xlabel('Longitud de onda (nm)')
    plt.ylabel('Intensidad')
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()

# Llamada a la función
plot_kind_class_outliners(class_outliners, X)

# imprime cuantos outliers de magnitud y forma hay
display(class_outliners['magnitude_outliers'])
print(f"Outliers de magnitud: {len(class_outliners['magnitude_outliers'])}")
print(f"Outliers de forma: {len(class_outliners['shape_outliers'])}")

array([ 10,  71, 131], dtype=int32)

Outliers de magnitud: 3
Outliers de forma: 22

Reproducir figura 2 del paper

Mediana funcional multivariada

Definición de la Profundidad Funcional Multivariada (MFD)

1. Proceso Considerado:

   • Se analiza un proceso estocástico $K$-variado $\{Y(t), t \in U\}$ definido en $\mathbb{R}^K$, con función de distribución acumulada $F_Y$ que genera caminos continuos en $C(U)^K$ (espacio de funciones continuas).

2. Medida de Profundidad:

   • La MFD para una curva arbitraria $X \in C(U)^K$ se define como: $$ \text{MFD}(X; F_Y) = \int_U D(X(t); F_{Y(t)}) \cdot w(t) \, dt, $$ donde:
   • $D$: Función de profundidad estadística en $\mathbb{R}^K$.
   • $w(t)$: Función de peso definida en $U$, que se normaliza a $1$, es decir, $\int_U w(t) dt = 1$.

3. Opciones para la Función de Peso:

Regiones de Profundidad

• Definición: La región de profundidad $ D_\alpha(F_X) $ en el nivel $ \alpha \geq 0 $ se define como: $$ D_\alpha(F_X) = \{ x \in \mathbb{R}^K : D(x; F_X) \geq \alpha \}, $$ donde $ D(x; F_X) $ es una función de profundidad que mide cuán central es $ x $ respecto a la distribución $ F_X $. Valores mayores de $ D(x; F_X) $ indican mayor centralidad.

Funciones de Peso

Las funciones de peso $ w(t) $ ajustan la contribución de diferentes puntos en el tiempo $ t $ en el análisis de datos funcionales. Se presentan dos ejemplos:

1. Peso Constante con Función Indicadora:

   • Un peso constante dentro de un rango de interés permite excluir ciertos periodos, como:
     • Fases de arranque en procesos industriales.
     • Mediciones imprecisas en regiones específicas de datos espectrales.

2. Peso Ajustado por Variabilidad:

   • Una función de peso que se adapta a la variabilidad local de los datos: $$ w(t) = w_\alpha(t; F_Y(t)) = \frac{\text{vol}(D_\alpha(F_Y(t)))}{\int_U\text{vol}(D_\alpha(F_{Y(u)})) \, du} $$ donde $ \text{vol}(D_\alpha(F_Y(t))) $ mide el "volumen" de la región de profundidad $ D_\alpha(F_Y(t)) $ en el tiempo $ t $. Este enfoque considera la variabilidad en amplitud (variabilidad vertical) de los datos funcionales.

---

Definición de Funciones y Correspondencias

1. Función $H$:

   • Se define como $H : U \times \mathbb{R}^K \to [0, 1]$: $$ H(t, x) = D(x, F_Y(t)), $$ donde $D(x, F_Y(t))$ es la profundidad del vector $x$ respecto a la distribución $F_Y(t)$.

2. Correspondencia $\gamma$:

   • Para cada $t \in U$, $\gamma(t)$ es el conjunto de $x \in \mathbb{R}^K$ en los que se evalúa la profundidad. Este conjunto es compacto y no vacío.

3. Funciones auxiliares:

   • $\Gamma: U \to \mathbb{R}$ se define como: $$ \Gamma(t) = (-\infty, \max_{x \in \gamma(t)} D(x; F_Y(t))]. $$
   • $G(t)$ denota el gráfico de $H(t, \cdot)$ restringido a $\gamma(t)$, y su clausura en $\mathbb{R}^K \times \mathbb{R} \cup \{-\infty\}$ se denota como $G(t)$.

Teorema 2: Propiedades del Máximo de Profundidad (Existencia de la mediana)

Suposiciones:

• $H$ es semicontinua superior.
• $\gamma$ es una correspondencia compacta y semicontinua superior.
• $\Gamma$ es semicontinua inferior.

Conclusiones:

1. La función $t \mapsto \max_{x \in \gamma(t)} D(x; F_Y(t))$ es continua.
2. El conjunto de puntos donde la profundidad máxima se alcanza: $$ \Pi(t) = \text{argmax}_{x \in \gamma(t)} D(x; F_Y(t)) $$ es compacto, no vacío y semicontinuo superior. Si $\Pi(t)$ es unívoco, es continuo.
3. Para una curva $\vartheta(t) \in \Pi(t)$ y $w(t)$ una función de peso, se tiene: $$ \text{MFD}(X; F_Y) = \int_U D(\vartheta(t); F_Y(t)) \cdot w(t) \, dt. $$
4. Si $\vartheta \in C(U)^K$ y maximiza la MFD, entonces $\vartheta(t) \in \Pi(t)$ en cada $t$.

La profundidad resultante se denomina profundidad de semiespacio funcional multivariante (MFHD).

Definición de Muestras Finitas

Para observaciones de curvas multivariantes $\{Y^1(t_j), \dots, Y^N(t_j)\}$ en tiempos discretos $t_1, \dots, t_T$:

1. Profundidad Muestral Multivariante Funcional (MFD): $$ \text{MFD}_N(X) = \sum_{j=1}^T D(X(t_j); F_Y(t_j), N) \cdot W_j, $$ donde $W_j$ es un peso que depende de $t_j$ y $w(t)$.

2. Consistencia: Se demuestra que $\text{MFD}_N$ converge a $\text{MFD}$ cuando tanto $N$ como $T$ tienden a infinito.

La media se define como la curva más profunda que para MFD por teorema 2 existe

lambdas_3 = data[:, :, 4:] 
display(lambdas_3.shape)
# Cambia la dimension de lambdas_3 intercambia por (571, 268, 3)
lambdas_3_new = np.moveaxis(lambdas_3, 0, 1)
display(lambdas_3_new.shape)
%R -i lambdas_3_new

(268, 571, 3)

(571, 268, 3)

%%R
library("mrfDepth")

median_mf <- mrfDepth::mfdmedian(lambdas_3_new)
median_mf

$MFDmedian
             [,1]        [,2]       [,3]
  [1,]   2.105005   0.8049057   1.587502
  [2,]   2.187033   0.8262345   1.696407
  [3,]   2.337496   0.8500005   1.888822
  [4,]   2.543371   0.8977931   2.194584
  [5,]   2.808741   0.9457164   2.610831
  [6,]   3.109565   1.0069587   3.140136
  [7,]   3.441962   1.0796899   3.748493
  [8,]   3.791654   1.1650708   4.530298
  [9,]   4.187731   1.2702830   5.472409
 [10,]   4.602064   1.4289436   6.637180
 [11,]   5.038211   1.5979974   7.921381
 [12,]   5.506924   1.8293113   9.483342
 [13,]   5.986744   2.0747602  11.140498
 [14,]   6.473741   2.3481041  13.186854
 [15,]   6.968772   2.6563346  15.320844
 [16,]   7.523273   2.9698299  17.749692
 [17,]   8.128189   3.3621962  20.502031
 [18,]   8.711750   3.8187841  23.694310
 [19,]   9.354988   4.2993876  26.936582
 [20,]   9.998370   4.7912929  30.372750
 [21,]  10.604316   5.4053845  34.377443
 [22,]  11.369136   6.0161060  38.780952
 [23,]  12.138011   6.7123730  43.627747
 [24,]  12.886487   7.4217204  48.598987
 [25,]  13.738171   8.2579607  53.902744
 [26,]  14.575708   9.1407035  60.028082
 [27,]  15.450214  10.0540269  66.105291
 [28,]  16.374814  11.0654319  72.785523
 [29,]  17.244564  11.9764134  78.900961
 [30,]  18.312495  13.2166005  86.890952
 [31,]  19.022477  14.1664380  92.618006
 [32,]  20.116644  15.3374546 100.484675
 [33,]  20.946326  16.3419502 107.213824
 [34,]  22.117082  17.6666607 115.869443
 [35,]  23.123445  18.8453012 121.897584
 [36,]  24.257791  20.1597170 131.805563
 [37,]  25.279525  21.3354233 137.616299
 [38,]  26.525818  22.7435449 147.237739
 [39,]  27.481334  24.0629645 153.810175
 [40,]  28.934868  25.5482467 163.883248
 [41,]  30.132987  27.1378440 170.173566
 [42,]  31.489360  28.5480889 179.655211
 [43,]  32.566663  29.6490323 187.907950
 [44,]  33.689635  31.0696740 195.343560
 [45,]  35.002795  32.5226736 204.341586
 [46,]  36.070027  33.8932282 211.257541
 [47,]  37.405291  35.3533386 219.636274
 [48,]  38.315319  36.6344593 225.285407
 [49,]  39.733432  38.0722759 231.668833
 [50,]  40.716470  39.4451926 237.818857
 [51,]  41.880403  40.8178573 243.546478
 [52,]  42.755958  42.0456595 249.612161
 [53,]  43.729239  42.9736918 254.974255
 [54,]  44.832729  44.1988506 260.858615
 [55,]  45.769214  45.2852136 265.970954
 [56,]  46.692657  46.4081007 271.443775
 [57,]  47.645147  47.5497123 274.399993
 [58,]  48.479773  48.5803963 279.004629
 [59,]  49.302679  49.5993938 282.368861
 [60,]  50.023688  50.5665519 284.815622
 [61,]  50.642496  51.4793485 287.497801
 [62,]  51.385963  52.3781814 289.308648
 [63,]  52.001433  53.0314963 291.130580
 [64,]  52.465882  53.7300048 292.176310
 [65,]  53.070219  54.4448569 293.609875
 [66,]  53.414099  55.1985523 295.837189
 [67,]  54.230735  56.1210323 296.862386
 [68,]  54.715249  56.9990554 298.768200
 [69,]  55.508267  57.8469431 299.003691
 [70,]  56.411754  58.4437657 300.650862
 [71,]  57.126491  59.2600650 300.777101
 [72,]  57.734251  60.0987221 300.711123
 [73,]  58.382562  60.8723601 301.878735
 [74,]  59.104979  61.6448273 301.763763
 [75,]  59.765328  62.3160232 301.530338
 [76,]  60.864386  63.0506014 301.295772
 [77,]  61.668376  63.7861976 301.288113
 [78,]  62.464861  64.5770600 300.372193
 [79,]  63.372372  65.2212285 299.885478
 [80,]  64.371142  65.8720298 299.164187
 [81,]  65.209276  66.8523518 297.668518
 [82,]  66.444708  67.5971799 296.686262
 [83,]  67.169503  67.9478712 294.661078
 [84,]  68.190111  68.8649447 294.570672
 [85,]  68.939222  69.4609354 291.548062
 [86,]  70.113131  70.1205653 290.914961
 [87,]  71.308014  70.9753590 291.068861
 [88,]  72.046522  71.6525559 289.301283
 [89,]  73.522340  72.4178235 289.256922
 [90,]  74.553856  73.4023682 288.420932
 [91,]  75.677896  74.3414597 288.487399
 [92,]  77.170862  75.0079096 287.501528
 [93,]  78.772393  76.0035190 285.077186
 [94,]  79.837337  76.5418750 285.063474
 [95,]  81.120451  77.3082635 283.418505
 [96,]  82.622196  78.1044900 283.894074
 [97,]  83.995226  78.9419488 281.893248
 [98,]  85.160245  79.6803612 280.786553
 [99,]  86.268862  80.5397237 279.399074
[100,]  87.302130  81.4258515 279.237465
[101,]  88.718670  82.0262374 277.435263
[102,]  90.015457  82.9256042 276.858727
[103,]  91.442713  83.8977394 275.422417
[104,]  92.697041  84.5633156 275.782874
[105,]  94.202897  85.4767722 273.763984
[106,]  95.606896  86.3002161 274.466289
[107,]  97.050653  87.0737461 273.813269
[108,]  98.432641  88.0843678 272.783111
[109,]  99.549037  88.7392417 271.490510
[110,] 100.975330  89.4130086 270.759227
[111,] 102.340477  90.1847089 269.369980
[112,] 103.658406  90.6050729 267.694011
[113,] 104.772167  91.1457487 266.702638
[114,] 106.193843  91.8482942 265.076138
[115,] 107.112255  92.3406410 263.160936
[116,] 108.356982  92.8690691 262.682393
[117,] 109.393185  93.5200331 260.829866
[118,] 110.582410  94.0708591 259.663137
[119,] 111.728098  94.6709009 259.688548
[120,] 112.860839  95.2643751 259.362386
[121,] 113.708150  95.4860430 258.483040
[122,] 115.293849  96.2248730 256.820611
[123,] 115.781633  96.5562556 256.593981
[124,] 116.458086  96.9907500 254.959935
[125,] 117.323280  97.4529138 253.843093
[126,] 118.531773  97.8266076 252.890896
[127,] 119.570154  98.3498613 252.299779
[128,] 120.367064  98.6225334 251.138060
[129,] 121.241266  99.1094153 250.870792
[130,] 121.835021  99.3150120 249.447604
[131,] 122.697827  99.5929174 249.070554
[132,] 122.907830  99.5181049 247.592968
[133,] 123.691463  99.9772003 247.127022
[134,] 124.531578 100.3340044 246.771631
[135,] 124.865036 100.4246486 245.558552
[136,] 125.334134 100.6581817 244.648907
[137,] 125.717981 100.6112761 243.923803
[138,] 126.270852 100.9430542 242.730814
[139,] 126.633448 101.1598014 242.061484
[140,] 126.976057 101.3887108 241.154035
[141,] 127.450503 101.5197168 240.537018
[142,] 127.950741 101.4732367 239.801896
[143,] 128.372640 101.8668673 239.261676
[144,] 128.682900 101.8953378 238.077881
[145,] 128.826146 101.8710395 237.458382
[146,] 128.740573 101.9278875 236.270508
[147,] 129.096224 102.1454971 235.586336
[148,] 129.305024 102.1706923 234.748144
[149,] 129.615068 102.3706690 234.024475
[150,] 129.829666 102.3768041 233.263806
[151,] 129.965066 102.4465602 232.525054
[152,] 130.028368 102.5625015 231.895422
[153,] 130.184601 102.6410795 231.059602
[154,] 130.318059 102.7343891 230.503040
[155,] 130.956052 103.0710442 230.224164
[156,] 130.954831 103.1024725 229.469580
[157,] 130.855169 103.1192938 228.568532
[158,] 130.622216 102.9142646 227.575026
[159,] 130.862398 103.2085299 227.260008
[160,] 130.835379 103.1736036 226.545650
[161,] 130.664783 102.9567262 225.752559
[162,] 130.312057 102.8380172 224.990103
[163,] 130.300527 102.8519813 224.516414
[164,] 130.521593 103.3841075 224.303310
[165,] 131.333597 104.0278120 224.287858
[166,] 130.585639 103.4720112 222.829166
[167,] 130.418296 103.6329203 222.137526
[168,] 129.413526 103.2753928 222.514421
[169,] 129.045711 103.1047861 221.136262
[170,] 128.726833 103.0796486 219.958045
[171,] 128.803166 103.1895023 218.911147
[172,] 128.233996 103.1147451 218.073913
[173,] 127.693202 103.0655587 217.130998
[174,] 127.417675 103.0134163 215.526439
[175,] 126.835367 102.6298819 215.494891
[176,] 125.881504 102.6107986 214.936164
[177,] 125.322834 102.5139776 214.045149
[178,] 125.002858 102.3793994 212.235979
[179,] 124.742566 102.2809006 210.825262
[180,] 124.409985 102.2239233 209.700114
[181,] 123.417046 102.1853923 209.162352
[182,] 122.674101 101.9789017 208.718170
[183,] 121.820939 101.8772595 207.963871
[184,] 121.724441 101.6641297 206.195380
[185,] 121.020380 101.5491043 205.097799
[186,] 120.910624 101.6130157 203.361406
[187,] 120.276186 101.2625467 202.059486
[188,] 119.370869 101.1271802 201.569717
[189,] 119.063494 100.9183789 200.547428
[190,] 118.167734 100.7582471 198.732096
[191,] 116.792695 100.5839388 198.949333
[192,] 116.497655 100.5694041 198.332284
[193,] 115.527783 100.6114467 196.932517
[194,] 115.024928 100.4023104 195.179109
[195,] 114.417714 100.2613160 194.731465
[196,] 113.909470 100.2736321 193.710491
[197,] 113.434083  99.9753197 192.559474
[198,] 112.821986  99.8383121 190.763013
[199,] 111.537973  99.7128760 189.657999
[200,] 110.289522  99.5556736 189.331370
[201,] 110.077287  99.6041187 187.443364
[202,] 108.887043  99.3398761 187.363985
[203,] 109.207405  99.3477490 184.903670
[204,] 108.664330  98.9232762 183.686359
[205,] 107.729279  98.8654745 183.048823
[206,] 106.711420  98.6172655 181.843206
[207,] 105.970114  98.6410924 181.056830
[208,] 105.335738  98.5524738 179.280944
[209,] 104.221403  98.3085533 179.213661
[210,] 103.280753  97.8924620 178.068998
[211,] 102.550153  97.8395338 177.181684
[212,] 102.614671  97.6923930 175.273280
[213,] 101.827870  97.6018089 173.724035
[214,] 100.901485  97.2632299 173.401112
[215,] 100.201867  97.0555162 172.076492
[216,]  99.352513  97.0163491 170.921400
[217,]  98.593180  96.7882731 170.159733
[218,]  97.853335  96.6105136 169.242697
[219,]  97.466315  96.6417549 167.919909
[220,]  96.548521  96.1799468 166.396024
[221,]  95.518471  96.1145910 166.423632
[222,]  94.845758  96.0468343 165.311850
[223,]  94.523744  96.1512785 164.596490
[224,]  93.513806  96.0175943 164.484434
[225,]  92.779019  95.7512629 163.371182
[226,]  92.540550  95.6970929 161.858189
[227,]  92.218382  95.5928066 160.886478
[228,]  91.267649  95.5047679 160.467959
[229,]  90.775148  95.4820795 159.705582
[230,]  90.565735  95.6951957 159.338035
[231,]  89.846931  95.5787396 158.736440
[232,]  89.293953  95.6129115 158.274820
[233,]  88.843813  95.4286769 157.656622
[234,]  88.234318  95.3500693 156.954268
[235,]  87.886669  95.5311066 156.386035
[236,]  87.388184  95.7231352 155.875580
[237,]  87.348655  95.7825825 154.837452
[238,]  86.506522  95.8437068 154.865861
[239,]  86.361072  95.7094117 153.811236
[240,]  85.519379  95.7326171 154.009034
[241,]  85.313088  95.6683310 153.041002
[242,]  84.719943  95.6682763 152.819855
[243,]  84.055191  95.4482427 151.989442
[244,]  83.845519  95.5570044 151.174997
[245,]  83.041307  95.6727176 150.734613
[246,]  82.763247  95.4405915 150.048871
[247,]  82.206473  95.0524500 149.271784
[248,]  81.778340  95.4621719 149.381010
[249,]  81.335401  95.3453500 148.520721
[250,]  81.403587  95.4250837 147.542226
[251,]  80.722578  95.3854522 146.858745
[252,]  80.833523  95.2775211 146.182450
[253,]  80.253653  95.3199663 145.829431
[254,]  79.740133  94.9237396 144.578935
[255,]  79.435585  95.0532062 144.103742
[256,]  78.829549  94.8869302 143.252742
[257,]  78.451354  94.6783548 142.614271
[258,]  77.761259  94.7945615 142.764083
[259,]  77.258806  94.5848047 142.244019
[260,]  77.100983  94.5875330 141.192032
[261,]  76.968030  94.6263475 140.247166
[262,]  76.521581  94.3024009 139.767875
[263,]  75.999684  94.5562247 139.506416
[264,]  75.891508  94.4147972 138.304081
[265,]  75.655621  94.6528889 138.097996
[266,]  75.594235  94.5072096 137.281704
[267,]  75.165432  94.4074324 136.909102
[268,]  74.970210  94.6429197 136.874606
[269,]  74.598544  94.4201499 136.571332
[270,]  74.419679  94.3898679 135.623790
[271,]  74.084147  94.5705727 135.430380
[272,]  73.647409  94.4719941 134.780073
[273,]  73.623847  94.4853537 134.274609
[274,]  73.475858  94.6956620 133.812457
[275,]  73.263236  94.7332009 133.537729
[276,]  73.204971  94.7638937 133.373090
[277,]  72.974472  94.8572602 133.110048
[278,]  72.998243  95.2735880 132.720986
[279,]  72.738265  95.2116584 132.906944
[280,]  72.744121  95.5105281 132.348176
[281,]  72.638573  95.6663775 131.808173
[282,]  72.823548  95.9463536 131.570332
[283,]  72.556854  95.9556729 131.147319
[284,]  72.612351  96.0811926 131.042762
[285,]  72.440545  96.1602553 130.613987
[286,]  72.279979  96.1612475 130.464944
[287,]  72.165040  96.4216479 130.226819
[288,]  71.776503  96.2299095 130.075279
[289,]  71.953473  96.5118170 129.537512
[290,]  71.764427  96.6551087 129.620354
[291,]  71.760978  96.8662763 129.308276
[292,]  71.798116  97.0492060 128.704293
[293,]  71.543434  96.8312613 128.046028
[294,]  71.346362  96.8805494 127.690261
[295,]  71.412662  96.9170203 127.201880
[296,]  71.290508  96.8338554 126.573302
[297,]  71.102466  96.7523767 126.265234
[298,]  70.805117  96.7847298 125.635177
[299,]  70.791130  96.7904069 125.444533
[300,]  70.868267  96.4732788 125.421771
[301,]  70.493065  96.3087908 125.112982
[302,]  70.404410  96.4964430 124.522419
[303,]  70.531671  96.7937439 124.128940
[304,]  70.602699  96.7109476 123.821249
[305,]  70.334441  96.8692391 123.550471
[306,]  70.170046  96.8043253 123.052652
[307,]  70.121549  96.7478605 122.755209
[308,]  70.354025  96.7601385 122.234825
[309,]  70.253305  97.0705875 122.160561
[310,]  70.135896  97.1464845 121.913085
[311,]  69.874392  97.0163407 121.258631
[312,]  69.785899  97.0696455 121.110859
[313,]  70.003710  97.1712884 120.755678
[314,]  69.848822  97.1963879 120.401093
[315,]  69.922598  97.4275375 120.170650
[316,]  69.927257  97.6765706 119.890643
[317,]  69.775766  97.8146695 119.601217
[318,]  69.865381  97.9265590 119.705000
[319,]  69.939432  98.1633412 119.520068
[320,]  70.013134  97.7196310 119.150946
[321,]  70.240042  97.8197485 118.964301
[322,]  70.088034  98.3800204 118.778117
[323,]  70.176711  97.9893529 118.486551
[324,]  70.607859  98.5796317 118.670520
[325,]  70.533625  98.5022533 118.355086
[326,]  70.545206  98.6435829 118.376397
[327,]  70.591806  98.6490788 118.082394
[328,]  70.668438  98.7125502 117.766900
[329,]  70.837114  98.8184388 117.385214
[330,]  70.690716  98.7074804 116.833833
[331,]  70.646488  98.8648709 117.119237
[332,]  70.590418  98.6900362 116.857496
[333,]  70.766485  98.8065552 116.848736
[334,]  70.805570  98.6694389 116.283979
[335,]  70.778672  98.7607651 116.180466
[336,]  70.976094  99.1230099 116.245108
[337,]  70.891753  98.9506982 115.762496
[338,]  70.902807  98.8635402 115.108349
[339,]  70.698457  98.6992138 115.014599
[340,]  70.846168  98.3744105 114.664121
[341,]  70.462085  98.2095985 113.933103
[342,]  70.283082  98.1509037 113.453540
[343,]  70.454126  98.2933285 112.691644
[344,]  70.243580  97.7580172 112.129044
[345,]  70.193647  97.6211087 111.958556
[346,]  70.078341  97.7549259 111.852585
[347,]  70.134841  97.3919161 111.198859
[348,]  69.591741  96.9765731 110.621779
[349,]  69.555315  96.6243819 109.772521
[350,]  69.150043  96.3762015 109.520727
[351,]  69.487061  96.2033717 108.788584
[352,]  68.794998  95.9360967 108.599812
[353,]  68.785639  95.9470727 108.250111
[354,]  68.438332  95.4908172 107.577582
[355,]  68.211680  95.2043990 107.127901
[356,]  67.812791  94.5537628 106.210793
[357,]  67.465485  94.1198868 105.561330
[358,]  67.363219  93.8823817 104.877522
[359,]  67.102348  93.1326387 104.287023
[360,]  66.958672  93.0052612 103.974080
[361,]  67.141048  92.8064252 103.062775
[362,]  66.918342  91.9980705 102.229498
[363,]  66.730119  91.9818140 101.532140
[364,]  66.463332  91.7456215 101.308491
[365,]  66.214728  91.0710338 100.535247
[366,]  65.912155  90.8775234 100.436073
[367,]  65.937461  90.5597000  99.340639
[368,]  65.835298  90.5662078  99.374654
[369,]  65.340586  90.1352563  98.787900
[370,]  65.116993  89.5522571  98.386044
[371,]  65.124216  89.4822619  98.030455
[372,]  64.560702  89.1954701  97.614758
[373,]  64.867071  88.6806681  96.962101
[374,]  64.483320  88.6105928  96.397552
[375,]  64.353096  88.5642197  96.158607
[376,]  64.214050  88.0972009  95.195176
[377,]  63.677964  87.4613327  94.554386
[378,]  63.304845  87.3430649  94.433003
[379,]  63.185703  87.1278510  94.155431
[380,]  63.072980  86.7685313  93.813315
[381,]  62.940424  86.5431471  93.337390
[382,]  62.613294  86.1487267  92.886841
[383,]  62.398085  85.7995405  92.325649
[384,]  62.324354  85.5257102  91.946834
[385,]  62.420648  85.2442565  91.204008
[386,]  62.054359  84.8125722  90.599992
[387,]  61.813623  84.4087462  90.217501
[388,]  61.136849  83.8625053  89.711102
[389,]  61.386586  83.4209744  88.998335
[390,]  60.924325  83.5447727  88.745019
[391,]  60.812772  82.7265793  88.161667
[392,]  60.290445  82.6078116  87.610872
[393,]  60.172295  81.8836750  86.656101
[394,]  59.550861  81.3429669  86.072844
[395,]  59.377899  81.0137760  85.530591
[396,]  59.307290  80.3465812  84.771780
[397,]  58.971026  80.2087495  84.550330
[398,]  58.992973  79.8403035  84.097634
[399,]  58.208248  78.9046891  83.049053
[400,]  58.111814  78.4739719  82.556720
[401,]  57.539510  78.1671024  82.082089
[402,]  57.225432  77.8345172  81.453342
[403,]  56.829731  77.2607806  80.811694
[404,]  56.316228  76.5923070  80.174146
[405,]  56.084056  75.9290886  79.573222
[406,]  55.500846  75.5724531  79.236260
[407,]  55.223778  74.9373075  78.617765
[408,]  54.947377  74.4694775  77.999512
[409,]  54.723625  73.9806539  77.361692
[410,]  54.336448  73.3371910  76.697368
[411,]  53.819250  72.9158745  75.905971
[412,]  53.477569  72.1481178  75.018336
[413,]  53.089046  71.8090646  74.676531
[414,]  52.787577  71.4121376  73.999247
[415,]  52.621226  71.0932569  73.513707
[416,]  51.996956  70.2681192  72.698755
[417,]  51.702334  69.7560198  72.171777
[418,]  51.263349  69.3132509  71.543478
[419,]  50.999568  68.6450940  70.531429
[420,]  50.647057  68.0530806  69.839189
[421,]  50.273899  67.6058682  69.483991
[422,]  49.742837  67.0734402  69.074112
[423,]  49.315085  66.5613991  68.593792
[424,]  48.906282  65.9481205  67.931353
[425,]  48.869655  65.4877590  67.318062
[426,]  48.408860  65.1133096  66.848963
[427,]  48.093356  64.4522535  66.122402
[428,]  47.885292  64.2250814  65.837219
[429,]  47.552243  63.6491926  65.430450
[430,]  47.198983  63.3085580  64.929072
[431,]  46.976763  62.8463555  63.866345
[432,]  46.506963  62.2155484  63.491783
[433,]  46.160509  61.6469987  63.099265
[434,]  45.625284  61.0133060  62.236673
[435,]  45.528323  60.7179898  62.167330
[436,]  45.075531  59.8908052  61.425378
[437,]  44.859211  59.6915714  61.103446
[438,]  44.522991  59.0415037  60.295042
[439,]  44.236413  58.8182428  59.970745
[440,]  43.902635  58.2198016  58.998117
[441,]  43.482127  57.8072981  58.759916
[442,]  43.008733  57.4327706  58.530083
[443,]  42.916562  57.0910084  58.065906
[444,]  42.736564  56.7135270  57.609575
[445,]  42.420203  56.3146550  56.929399
[446,]  42.197105  55.8569975  56.610375
[447,]  41.880534  55.1243086  55.915470
[448,]  41.439165  54.5437864  55.290622
[449,]  41.281269  54.5114432  55.124416
[450,]  40.848898  53.7409114  54.453948
[451,]  40.407669  53.2586410  54.092018
[452,]  40.260678  52.8520817  53.456213
[453,]  39.745399  52.4605903  52.927159
[454,]  39.437811  52.0027810  52.485045
[455,]  39.063678  51.4410683  52.116032
[456,]  38.845305  50.9189721  51.549922
[457,]  38.619696  50.6058922  51.163774
[458,]  38.408089  50.3478571  50.595062
[459,]  38.138368  49.7591114  50.177244
[460,]  37.636884  49.2559570  49.615593
[461,]  37.146041  48.7490251  49.307074
[462,]  36.843877  48.2950388  48.620835
[463,]  36.484139  47.9623060  47.992868
[464,]  36.132799  47.2826583  47.597230
[465,]  35.822218  46.7906017  47.039541
[466,]  35.424274  46.2019437  46.449483
[467,]  35.366665  45.7672959  45.914468
[468,]  34.944337  45.2994888  45.488368
[469,]  34.557196  44.8253570  45.118191
[470,]  34.181615  44.3431641  44.538207
[471,]  33.997883  43.9219635  44.130198
[472,]  33.527514  43.4300813  43.642760
[473,]  33.037279  42.5837760  42.785981
[474,]  32.857912  42.5492558  42.624348
[475,]  32.338975  41.6157093  41.989899
[476,]  32.298649  41.2597636  41.527587
[477,]  31.915716  40.9048877  41.068086
[478,]  31.649197  40.4520438  40.549889
[479,]  31.346764  39.7309124  39.956620
[480,]  30.994531  39.3971521  39.566051
[481,]  30.995270  38.9641738  39.140716
[482,]  30.650329  38.5899307  38.780284
[483,]  30.392964  38.0516385  38.200803
[484,]  30.368529  37.6689670  37.796375
[485,]  30.150901  37.1393331  37.181001
[486,]  29.850110  36.6162946  36.739157
[487,]  29.673507  36.1295034  36.352445
[488,]  29.512182  35.8957968  35.946818
[489,]  29.377477  35.3518461  35.499205
[490,]  29.395687  34.9292181  35.122826
[491,]  29.091301  34.4748099  34.654455
[492,]  29.007689  34.1722687  34.300845
[493,]  28.770854  33.8364482  33.907873
[494,]  28.464079  33.4699600  33.417387
[495,]  28.250202  33.0658648  33.047173
[496,]  28.031538  32.7155852  32.668705
[497,]  27.619892  32.2154810  32.166836
[498,]  27.507217  31.8911028  31.830635
[499,]  27.302721  31.6402086  31.579664
[500,]  27.004787  31.0184236  30.961136
[501,]  26.802900  30.8283357  30.754387
[502,]  26.426886  30.3641909  30.414630
[503,]  26.169907  30.2318988  30.220725
[504,]  25.826648  29.8890106  29.873011
[505,]  25.469211  29.6519594  29.595173
[506,]  25.298590  29.3125703  29.298179
[507,]  24.873730  28.9727597  28.912925
[508,]  24.416003  28.7932364  28.655387
[509,]  24.140214  28.5998858  28.326844
[510,]  23.799464  28.3522355  28.048803
[511,]  23.373861  28.1424942  27.794277
[512,]  22.997167  27.8283860  27.437554
[513,]  22.611955  27.5448211  27.159971
[514,]  22.342379  27.2180932  26.847463
[515,]  22.086065  27.0243366  26.601820
[516,]  21.726938  26.6793888  26.238083
[517,]  21.389398  26.4368822  25.980092
[518,]  21.145042  26.1512298  25.699041
[519,]  20.844910  25.9067654  25.464228
[520,]  20.490598  25.7985359  25.190087
[521,]  20.211713  25.5626125  24.970792
[522,]  20.005238  25.2971877  24.696478
[523,]  19.762824  25.1726801  24.510781
[524,]  19.499280  24.8999935  24.273504
[525,]  19.298056  24.7587017  24.039962
[526,]  19.034981  24.4386384  23.651369
[527,]  18.848639  24.2259704  23.433126
[528,]  18.721624  24.0030799  23.219574
[529,]  18.444675  23.7220301  22.903975
[530,]  18.255543  23.4595183  22.599627
[531,]  18.034231  23.1664772  22.302126
[532,]  17.979610  22.8303775  22.102634
[533,]  17.862578  22.5210602  21.923181
[534,]  17.617786  22.2619123  21.621950
[535,]  17.505179  22.0778776  21.416549
[536,]  17.349889  21.8512729  21.235292
[537,]  17.235709  21.5835744  21.014166
[538,]  17.122858  21.4647233  20.813460
[539,]  16.955003  21.1913191  20.492590
[540,]  16.824548  20.9156469  20.283671
[541,]  16.696310  20.6032265  19.993806
[542,]  16.599919  20.3462238  19.774708
[543,]  16.510182  20.0923079  19.567411
[544,]  16.354441  19.7008154  19.317520
[545,]  16.259869  19.4876605  19.138188
[546,]  16.228117  19.2508215  19.018714
[547,]  16.133068  19.0151594  18.830690
[548,]  16.036073  18.7421263  18.643073
[549,]  15.946948  18.5261823  18.432642
[550,]  15.877955  18.2749890  18.213944
[551,]  15.803750  17.9681025  17.979374
[552,]  15.730516  17.7996062  17.792235
[553,]  15.722362  17.5850027  17.628365
[554,]  15.634058  17.3849822  17.426879
[555,]  15.642303  17.1934073  17.309121
[556,]  15.594689  16.8640576  17.118360
[557,]  15.520770  16.7454197  16.960399
[558,]  15.477818  16.4796340  16.802293
[559,]  15.455418  16.3392477  16.733105
[560,]  15.433770  16.1771071  16.604196
[561,]  15.388948  15.9217118  16.458106
[562,]  15.419322  15.7590117  16.375666
[563,]  15.369182  15.5948802  16.289092
[564,]  15.376938  15.4044170  16.209636
[565,]  15.363430  15.2604665  16.139103
[566,]  15.360895  15.1452977  16.133887
[567,]  15.321499  14.9911914  16.140242
[568,]  15.354468  14.8634264  16.205643
[569,]  15.352377  14.7269768  16.260211
[570,]  15.422172  14.6197279  16.368587
[571,]  15.459469  14.4675093  16.401015

attr(,"class")
[1] "mrfDepth"  "mfdmedian"

Cargando paquete requerido: ggplot2

%R -o median_mf

import numpy as np
import matplotlib.pyplot as plt
from collections import OrderedDict

def plot_lambdas(lambdas_3, median_mf, n_samples=None):
    """
    Plot intensity data for lambdas_3 with median MFD curves
    
    Parameters:
    - lambdas_3: numpy array with shape (n_samples, points, wavelengths)
    - median_mf: OrderedDict with median MFD curves
    - n_samples: number of samples to plot (default: all samples)
    """
    if n_samples is None:
        n_samples = lambdas_3.shape[0]
    
    # Limit to specified number of samples
    lambdas_subset = lambdas_3[:n_samples]
    
    # Create color gray palette
    colors = plt.cm.gray(np.linspace(0, 1, n_samples))
    
    plt.figure(figsize=(15, 6))
    
    # Plot each wavelength for samples
    for wavelength in range(lambdas_subset.shape[2]):
        for sample in range(n_samples):
            plt.plot(lambdas_subset[sample, :, wavelength],
                     color=colors[sample],
                     linewidth=.5,
                     alpha=0.3)
    
    # Color palette for median curves
    labels = [r"Mediana de $\lambda = 5$" , r"Mediana de $\lambda = 6$", r"Mediana de $\lambda = 7$"]
    
    # Plot median MFD curves
    key, curve = next(iter(median_mf.items()))
    plt.plot(curve, linewidth=3, label= labels)
    plt.title(f'Lambdas Data - First {n_samples} Samples', fontsize=16)
    plt.xlabel('Measurement Points')
    plt.ylabel('Intensity')
    plt.grid(True, alpha=0.3)
    plt.legend()
    
    plt.tight_layout()
    plt.show()
# Ejemplo de uso (descomentar y ajustar según tus datos)
# plot_lambdas(lambdas_3, median_mf=median_mf)

# Ejemplo de uso 
plot_lambdas(lambdas_3, median_mf=median_mf)

Outliers del proceso funcional multivariado

Teniendo un vector de funciones $\mathbf{X}$

el primer paso es la profundidad como, tomando pointwise

$$\text{MFMD} = \int_T D(\mathbf{X} (t), F_{\mathbf{X} (t)})\omega(t) dt$$

Donde $D$ es una profundidad multivariada y MFMD es una profundidad funcional multivariada.

$$ D: \mathcal{R}^p \to \mathcal{R}$$ $$\text{MFMD}: \mathcal{L}^{2p}_T \to \mathcal{R}_{[0, 1]}$$

Se define una función de outlyingness como $$ o: \mathcal{R}^p \to [0 , \infty)$$

La profundidad funcional multivariada se puede definir con base a la función de outlyingness (atipicidad) como

$$ D(\mathbf{X} (t), F_{\mathbf{X} (t)}) = [1 + o(\mathbf{X} (t), F_{\mathbf{X} (t)})]^{-1}$$ Despejando:

$$ o(\mathbf{X} (t), F_{\mathbf{X} (t)}) = [D(\mathbf{X} (t), F_{\mathbf{X} (t)})^{-1} - 1] $$

Ahora se define la función O por medio de la O:

$$ O((\mathbf{X} (t), F_{\mathbf{X} (t)})) = o((\mathbf{X} (t), F_{\mathbf{X} (t)}))\vec{v}((\mathbf{X} (t), F_{\mathbf{X} (t)}))$$

$$ O: \mathcal{R}^p \to \mathcal{R}^p$$

Donde $\vec{v}(\mathbf{y})$ es un vector unitario en la dirección de $\mathbf{y}$

$$ \vec{v}(\mathbf{y}) = \frac{\mathbf{y} - \tilde{\mathbf{y}}}{||\mathbf{y} - \tilde{\mathbf{y}}||} = \text{Signo del vector y} = Sign(\mathbf{y})$$

Anomalía Direccional Funcional

Se define un proceso estocástico $ X : I \to \mathbb{R}^p $ que toma valores en un espacio de funciones continuas $ C(I, \mathbb{R}^p) $, con distribución de probabilidad $ F_X $. Para medir la anomalía direccional de $ X $, se utilizan las siguientes métricas:

1. Anomalía direccional funcional ($ FO $): Representa la anomalía total de $ X $, acumulando la magnitud de la anomalía direccional a lo largo del intervalo $ I $: $$ FO(X, F_X) = \int_I \|O(X(t), F_{X(t)})\|^2 w(t) \, dt $$ Aquí, $ O(X(t), F_{X(t)}) $ es la anomalía direccional en el tiempo $ t $, y $ w(t) $ es una función de peso.

2. Anomalía direccional media ($ MO $): Evalúa la posición promedio de $ X $ en relación con la curva central, considerando tanto la distancia como la dirección: $$ MO(X, F_X) = \int_I O(X(t), F_{X(t)}) w(t) \, dt $$ $$ =\int_I\left[ f , f , \dots, f \right]$$ $$ =\left[ \int_I , \int_I , \dots, \int_I \right] \in \mathcal{R}^p$$

• Note que MO devuelve un vector, (se esta integrando un vector), las funciones f pertenecen a $\mathcal{R}^p$
• La norma de MO da la la magnitud de la atipicidad que son outliers de magnitud.

3. Variación de la anomalía direccional ($ VO $): Mide las variaciones de $ O(X(t), F_{X(t)}) $ en términos de norma y dirección a lo largo del intervalo: $$ VO(X, F_X) = \int_I \|O(X(t), F_{X(t)}) - MO(X, F_X)\|^2 w(t) \, dt $$

• Note que VO es un número.
• Determina los outliners de forma.

Se puede definir la distancia de mahalanobis robusta como

$$ \vec{y} = (\underbrace{\vec{MO}}_{\mathcal{R}^p}, \underbrace{VO}_{\mathcal{R}}) \in \mathcal{R}^{p +1}$$

La función de covarianza cambia como se usa MCD (minimum covariance determinant) para encontrar la covarianza robusta.

$$ \frac{1}{h} \sum_{i=1}^h (\vec{y}_i - \bar{\vec{y}})(\vec{y}_i - \bar{\vec{y}})^T$$

Distribución y Detección de Valores Atípicos en Curvas Generadas por Procesos Gaussianos

A través de estudios numéricos, se encontró que la distribución de $ Y_{k,n} = (MO^T_{T_k,n}, VO_{T_k,n})^T $ se puede aproximar bien mediante una distribución normal de dimensión $ p+1 $ cuando $ X $ es generado por un proceso estacionario gaussiano de $ p $ dimensiones.

Detección de Valores Atípicos:

• Se utilizó la distancia de Mahalanobis robusta ($ RMD^2 $) junto con los estimadores mínimos del determinante de covarianza (MCD)

Procedimiento para Detectar Valores Atípicos

1. Cálculo de la distancia de Mahalanobis robusta: $$ RMD^2 (Y_{k,n}, \bar{Y}^*_{k,n,J}) = (Y_{k,n} - \bar{Y}^*_{k,n,J})^T S^*_{k,n,J}{}^{-1}(Y_{k,n} - \bar{Y}^*_{k,n,J}), $$ donde:

   • $ J $: grupo de $ h $ puntos que minimizan el determinante de la matriz de covarianza.
   • $ \bar{Y}^*_{k,n,J} = h^{-1} \sum_{i \in J} Y_{k,n,i} $.
   • $ S^*_{k,n,J} = h^{-1} \sum_{i \in J} (Y_{k,n,i} - \bar{Y}^*_{k,n,J})(Y_{k,n,i} - \bar{Y}^*_{k,n,J})^T $.
   • $ h $: submuestra que controla la robustez del método.

2. Aproximación de la cola de la distribución: Según Hardin y Rocke (2005), la cola de $ RMD^2 $ sigue: $$ \frac{c(m - p)}{m(p + 1)} RMD^2 (Y_{k,n}, \bar{Y}^*_{k,n,J}) \sim F_{p+1,m-p}, $$ donde:

   • $ c $ y $ m $: parámetros para la distribución $ F $ y el factor de escala, calculados mediante simulación.
   • El valor de corte ($ C $) es el cuantil $ \alpha $ ($ \alpha = 0.993 $) de $ F_{p+1,m-p} $.

3. Identificación de valores atípicos: Una curva se clasifica como atípica si: $$ \frac{c(m - p)}{m(p + 1)} RMD^2 (Y_{k,n}, \bar{Y}^*_{k,n,J}) > C. $$

Usos Adicionales del $ RMD $

• Medida de centralidad para definir:
  • La mediana y la región central de las curvas.
  • La función de media recortada.
  • El boxplot funcional.

lambdas_3.shape

(268, 571, 3)

%R -i lambdas_3

%%R 

library("fdaoutlier")
# Show mr the dimensions of lambdas_3

print(dim(lambdas_3))

directional_outliners <- dir_out(
  lambdas_3,
  data_depth = "random_projections",
  n_projections = 200L,
  seed = NULL,
  return_distance = TRUE,
  return_dir_matrix = TRUE
)

[1] 268 571   3

Con las distancias es necesario hacer re muestreo de n, con 1000 veces ajusto una distribución y los que esten fuera del quantil 99.3

%R -o directional_outliners
# show me the components of the directional_outliners 
display(directional_outliners.keys())
display(directional_outliners['mean_outlyingness'].shape)
display(directional_outliners['var_outlyingness'].shape)
display(directional_outliners['ms_matrix'].shape)
display(directional_outliners['distance'].shape)

(np.str_('mean_outlyingness'),
 np.str_('var_outlyingness'),
 np.str_('distance'),
 np.str_('ms_matrix'),
 np.str_('mcd_obj'),
 np.str_('dirout_matrix'))

(268, 3)

(268,)

(268, 4)

(268,)

distance = directional_outliners['distance']

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

def bootstrap_percentile_distribution(data, n_resamples=1000, percentile=99.3):
    """
    Realiza remuestreo bootstrap, grafica la distribución empírica 
    y marca el percentil especificado.
    
    Parámetros:
    -----------
    data : list or numpy.ndarray
        Lista o array de datos originales
    n_resamples : int, opcional (por defecto 1000)
        Número de remuestreos bootstrap
    percentile : float, opcional (por defecto 99.7)
        Percentil para determinar valores extremos
    
    Retorna:
    --------
    tuple
        - Array de índices de valores originales que superan el percentil
        - Percentil calculado
    """
    # Convertir los datos a un array de NumPy
    data = np.array(data)
    
    # Número de muestras
    n = len(data)
    
    # Almacenar resultados de los remuestreos
    bootstrap_samples = np.zeros((n_resamples, n))
    
    # Realizar remuestreo bootstrap
    for i in range(n_resamples):
        # Remuestreo con reemplazo
        bootstrap_samples[i] = np.random.choice(data, size=n, replace=True)
    
    # Calcular el percentil para cada remuestreo
    bootstrap_percentiles = np.percentile(bootstrap_samples, percentile, axis=1)
    
    # Calcular el cuantil de corte
    quantil_corte = np.percentile(bootstrap_percentiles, percentile)
    
    # Encontrar los índices de los valores originales por encima del percentil
    extreme_indices = np.where(data > quantil_corte)[0]
    
    # Configurar el estilo de la gráfica
    plt.figure(figsize=(10, 6))
    sns.set_style('whitegrid')
    
    # Graficar la distribución empírica de los remuestreos
    sns.histplot(bootstrap_percentiles, kde=True, color='blue', alpha=0.7)
    
    # Añadir línea vertical roja punteada para el cuantil de corte
    plt.axvline(x=quantil_corte, color='red', linestyle='--', 
                label=f'Cuantil {percentile}%: {quantil_corte:.4f}')
    
    plt.title('Distribución Empírica de Percentiles Bootstrap')
    plt.xlabel('Distancia robusta de Mahalanobis')
    plt.ylabel('Frecuencia')
    plt.legend()
    
    plt.tight_layout()
    plt.show()
    
    return extreme_indices, quantil_corte


indices_extremos, cuantil = bootstrap_percentile_distribution(distance)
    
print("Índices de valores extremos:", indices_extremos)
print(f"Cuantil de corte (99.3): {cuantil}")
print("Valores extremos:", distance[indices_extremos])

Índices de valores extremos: []
Cuantil de corte (99.3): 2443.7431231536584
Valores extremos: []

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

def estimate_distribution_parameters(distances, 
                                     n_bootstrap: int = 1000, 
                                     alpha: float = 0.993, 
                                     seed: int = None):
    """
    Estimate m and c parameters for the F-distribution based approach
    
    Parameters:
    -----------
    distances : np.ndarray or list
        Input distances to analyze
    n_bootstrap : int, optional
        Number of bootstrap iterations
    alpha : float, optional
        Significance level for outlier cutoff
    seed : int, optional
        Random seed for reproducibility
    
    Returns:
    --------
    Tuple of (estimated m, estimated c, standard deviation of c)
    """
    # Convert to numpy array if not already
    distances = np.asarray(distances)
    flat_distances = distances.flatten()
    
    if seed is not None:
        np.random.seed(seed)
    
    # Number of features (p)
    p = 1  # Assuming univariate case
    
    # Prepare arrays to store bootstrap results
    m_estimates = np.zeros(n_bootstrap)
    c_estimates = np.zeros(n_bootstrap)
    
    # Bootstrapping loop
    for i in range(n_bootstrap):
        # Resample distances with replacement
        bootstrap_sample = np.random.choice(flat_distances, size=len(flat_distances), replace=True)
        
        # Try different m values
        m_candidates = np.linspace(len(flat_distances)/2, len(flat_distances)*2, 20)
        best_error = np.inf
        best_m = len(flat_distances)
        best_c = 1.0
        
        for m_candidate in m_candidates:
            # Scale distances according to the formula
            scaled_distances = bootstrap_sample * (m_candidate - p) / (m_candidate * (p + 1))
            
            try:
                # Theoretical F-distribution quantile
                theoretical_quantile = stats.f.ppf(alpha, p + 1, m_candidate - p)
                
                # Empirical quantile
                empirical_quantile = np.quantile(scaled_distances, alpha)
                
                # Calculate error
                error = abs(theoretical_quantile - empirical_quantile)
                
                # Update best parameters if error is smaller
                if error < best_error:
                    best_error = error
                    best_m = m_candidate
                    best_c = theoretical_quantile / empirical_quantile
            except Exception:
                continue
        
        # Store the best estimates
        m_estimates[i] = best_m
        c_estimates[i] = best_c
    
    # Remove any zero or invalid estimates
    m_estimates = m_estimates[m_estimates > 0]
    c_estimates = c_estimates[c_estimates > 0]
    
    # Compute statistics
    mean_m = np.mean(m_estimates)
    mean_c = np.mean(c_estimates)
    std_c = np.std(c_estimates)
    
    return mean_m, mean_c, std_c

def detect_outliers(distances, 
                    alpha: float = 0.993, 
                    bootstrap_m_c: bool = True,
                    n_bootstrap: int = 1000,
                    plot: bool = True):
    """
    Detect outliers using robust distance method
    
    Parameters:
    -----------
    distances : np.ndarray or list
        Input distances to analyze
    alpha : float, optional
        Significance level for outlier cutoff
    bootstrap_m_c : bool, optional
        Whether to estimate m and c via bootstrapping
    n_bootstrap : int, optional
        Number of bootstrap iterations
    plot : bool, optional
        Whether to generate diagnostic plots
    
    Returns:
    --------
    List of outlier indices
    """
    # Convert to numpy array if not already
    distances = np.asarray(distances)
    flat_distances = distances.flatten()
    
    # Number of features (p)
    p = 1  # Assuming univariate case
    
    # Estimate m and c via bootstrapping if requested
    if bootstrap_m_c:
        m, c_factor, c_std = estimate_distribution_parameters(
            distances, 
            n_bootstrap=n_bootstrap, 
            alpha=alpha
        )
    else:
        # Default values if not bootstrapping
        m = len(flat_distances)
        c_factor = 1.0
        c_std = None
    
    # Scale robust distances
    scaled_distances = c_factor * (m - p) / (m * (p + 1)) * flat_distances
    
    # Cutoff value from F-distribution
    cutoff_value = stats.f.ppf(alpha, p + 1, m - p)
    
    # Detect outliers
    is_outlier = scaled_distances > cutoff_value
    
    # Create diagnostic plots if requested
    if plot:
        plt.figure(figsize=(15, 5))
        
        # Original Distances
        plt.subplot(1, 3, 1)
        plt.hist(flat_distances, bins=30, density=True, alpha=0.7)
        plt.title('Original Distances')
        plt.xlabel('Distance')
        
        # Scaled Distances
        plt.subplot(1, 3, 2)
        plt.hist(scaled_distances, bins=30, density=True, alpha=0.7)
        plt.title('Scaled Distances')
        plt.axvline(cutoff_value, color='r', linestyle='--', label='Cutoff')
        plt.xlabel('Scaled Distance')
        plt.legend()
        
        # Box Plot of Distances
        plt.subplot(1, 3, 3)
        plt.boxplot(flat_distances)
        plt.title('Boxplot of Distances')
        plt.xlabel('Distance')
        
        plt.tight_layout()
        plt.show()
        
        # Q-Q Plot
        plt.figure(figsize=(10, 5))
        stats.probplot(flat_distances, dist="f", sparams=(p+1, m-p), plot=plt)
        plt.title('Q-Q Plot of Distances')
        plt.show()
        
        # Print diagnostic information
        print(f"Estimated m: {m:.4f}")
        print(f"Estimated c factor: {c_factor:.4f}")
        if c_std is not None:
            print(f"Standard deviation of c: {c_std:.4f}")
        print(f"Cutoff value: {cutoff_value:.4f}")
        print(f"Number of outliers: {sum(is_outlier)}")
    
    # Return indices of outliers
    return np.where(is_outlier)[0].tolist()

print("\n--- Without Bootstrapping M and C ---")
outliers_without_bootstrap = detect_outliers(distance, bootstrap_m_c=False)

print("\n--- With Bootstrapping M and C ---")
outliers_with_bootstrap = detect_outliers(distance, bootstrap_m_c=True)
print(f"Los outliers detectados son: {outliers_with_bootstrap}")
# con retorno

--- Without Bootstrapping M and C ---

Estimated m: 268.0000
Estimated c factor: 1.0000
Cutoff value: 5.0552
Number of outliers: 103

--- With Bootstrapping M and C ---

Estimated m: 134.0000
Estimated c factor: 0.0118
Standard deviation of c: 0.0050
Cutoff value: 5.1516
Number of outliers: 3
Los outliers detectados son: [70, 130, 267]

Tareas Taller 1

Profundidad de variación total y similitud de forma

import os
import numpy as np
import skfda
from skfda.representation.basis import BSplineBasis
from skfda.preprocessing.smoothing import BasisSmoother
import matplotlib.pyplot as plt
from tqdm import tqdm

def smooth_functional_data(
    data, 
    wavelengths, 
    dominio, 
    n_basis=65, 
    best_param=1e-3, 
    rango_curvas=[None, None]
):
    # Determinar el rango de curvas a procesar
    inicio = rango_curvas[0] if rango_curvas[0] is not None else 0
    fin = rango_curvas[1] if rango_curvas[1] is not None else data.shape[0]

    # Inicializar array para almacenar datos suavizados
    smoothed_data = np.zeros_like(data[inicio:fin])

    # Crear una barra de progreso
    with tqdm(total=(fin-inicio)*len(wavelengths), desc="Suavizando muestras") as pbar:
        # Para cada longitud de onda
        for wave_idx, _ in wavelengths.items():
            # Procesar muestras para esta longitud de onda
            for i, muestra in enumerate(range(inicio, fin)):
                # Obtener los datos para esta muestra y longitud de onda
                X = data[muestra, :, wave_idx].flatten()
                
                # Crear una base B-Spline
                basis = BSplineBasis(domain_range=(275, 560), n_basis=n_basis)
                
                # Crear functional data grid
                fd = skfda.FDataGrid(data_matrix=X, grid_points=dominio)
                
                # Suavizado penalizado
                smoother = BasisSmoother(basis, smoothing_parameter=best_param)
                fd_smoothed = smoother.fit_transform(fd)
                
                # Guardar los datos suavizados (flatten to ensure correct shape)
                smoothed_data[i, :, wave_idx] = fd_smoothed(dominio)[0].flatten()
                
                # Actualizar la barra de progreso
                pbar.update(1)

    return smoothed_data

init = 5*268
end = 6*268

# Ejemplo de uso
smoothed_result = smooth_functional_data(data, wavelengths, dominio)

Suavizando muestras: 100%|██████████| 1876/1876 [05:30<00:00,  5.68it/s]

print(class_outliners.keys())

# Obtener los outliers de forma y magnitud
tvd = class_outliners['tvd']
mss = list(class_outliners['mss'])
shp_out = list(class_outliners['shape_outliers'])
# has que shp_out no tenga componentes del tipo np.int32
shp_out = [int(i) for i in shp_out]
mg_out = list(class_outliners['magnitude_outliers'])

(np.str_('outliers'), np.str_('shape_outliers'), np.str_('magnitude_outliers'), np.str_('tvd'), np.str_('mss'))

1. Boxplot de MSS

import matplotlib.pyplot as plt
import numpy as np

def boxplot_mss(data, f=1.5):
    """
    Crea un boxplot y detecta los outliers de una lista de valores.
    
    Args:
        data (list): Lista de valores numéricos
        f (float): Factor para calcular los límites de los outliers (default: 1.5)
        
    Returns:
        list: Lista con los índices de los outliers encontrados
    """
    # Crear la figura
    plt.figure(figsize=(8, 6))
    
    # Crear el boxplot
    bp = plt.boxplot(data, patch_artist=False)
    
    # Calcular estadísticas
    Q1 = np.percentile(data, 25)
    Q3 = np.percentile(data, 75)
    IQR = Q3 - Q1
    lower_bound = Q1 - f * IQR
    upper_bound = Q3 + f * IQR
    
    # Identificar índices de outliers
    outlier_indices = [i for i, x in enumerate(data) if x < lower_bound or x > upper_bound]
    outlier_values = [data[i] for i in outlier_indices]
    
    # Plotear los outliers en color naranja
    if len(outlier_values) > 0:
        plt.plot([1] * len(outlier_values), outlier_values, 'o', color='orange', label='Outliers')
    
    # Configurar el gráfico
    plt.title('Shape outlyingness Boxplot')
    plt.grid(True, linestyle='--', alpha=0.7)
    
    # Añadir texto con los límites
    plt.figtext(0.02, 0.02, 
                f'Límite inferior: {lower_bound:.2f}\n'
                f'Límite superior: {upper_bound:.2f}', 
                fontsize=8)
    
    # Mostrar el gráfico
    plt.show()
    
    return outlier_indices

shape_out = boxplot_mss(mss, 3)
# a todos los valores de shape_out súmalele 1
shape_out_correction = [i+1 for i in shape_out]
print(f"Índices de los outliers (py): {shape_out_correction}")

print(f"Índices de los outliers (R ): {shp_out}")

Índices de los outliers (py): [104, 109, 125, 129, 138, 140, 143, 145, 147, 149, 150, 151, 152, 153, 154, 155, 159, 160, 162, 163, 164, 167]
Índices de los outliers (R ): [104, 109, 125, 129, 138, 140, 143, 145, 147, 149, 150, 151, 152, 153, 154, 155, 159, 160, 162, 163, 164, 167]

2. functional boxplot TVD

def functional_boxplot_tvd(X, shape_out, F=1.5):
    """
    Create functional boxplot with outliers
    """
    # Filtra las filas tal que NO esten en shape_out
    X_out = X
    fd = skfda.FDataGrid(data_matrix=X_out, grid_points=dominio, dataset_name=f'Functional Boxplot TVD')
    fdBoxplot = Boxplot(fd, factor=F)        
    fdBoxplot.plot()
    return fdBoxplot
    

fboxplot = functional_boxplot_tvd(X, shape_out, 1.5)
outliers = fboxplot.outliers
# get the index where the outliers are
index_out = list(np.where(outliers)[0])

magnitude_out = [int(x)  for x in index_out]

index_out = [x + 1 for x in magnitude_out]

mg_out = [int(x) for x in mg_out]

print(f"Índices de los outliers (py): {index_out}")
print(f"Índices de los outliers (R ): {mg_out}")

Índices de los outliers (py): [10, 71, 131]
Índices de los outliers (R ): [10, 71, 131]

import numpy as np
import matplotlib.pyplot as plt

def plot_kind_class_outliners(magnitude_outliers, shape_outliers, X):
    """
    Plot outliers and non-outliers with different colors and labels.
    Shows only 3 shape outliers and all other curves in grey with alpha=0.5
    
    Args:
        magnitude_outliers (list): List of indices for magnitude outliers
        shape_outliers (list): List of indices for shape outliers
        X (numpy.ndarray): Array with all curves data
    """
    # Convierte los índices de outliers a conjuntos para una búsqueda eficiente
    magnitude_outliers_set = set(magnitude_outliers)
    shape_outliers_set = set(shape_outliers)
    
    # Crear máscaras para diferentes tipos de outliers
    magnitude_outliers_mask = np.array([i in magnitude_outliers_set for i in range(len(X))])
    shape_outliers_mask = np.array([i in shape_outliers_set for i in range(len(X))])
    
    # Encuentra los magnitude outliers en X
    magnitude_curves = X[magnitude_outliers_mask, :]
    
    # Encuentra los shape outliers en X (limitado a 3)
    shape_curves = X[shape_outliers_mask, :]
    shape_curves = shape_curves[:3] if len(shape_curves) > 3 else shape_curves
    
    # Encuentra las curvas que no son outliers
    non_outliers_mask = ~(magnitude_outliers_mask | shape_outliers_mask)
    non_outliers = X[non_outliers_mask, :]
    
    # Grafica los tipos de curvas
    plt.figure(figsize=(12, 6))
    
    # Primero grafica todas las curvas no outliers en gris
    for i in range(non_outliers.shape[0]):
        plt.plot(non_outliers[i], color='gray', alpha=0.5, label='_')
    
    # Grafica los magnitude outliers
    for i in range(magnitude_curves.shape[0]):
        plt.plot(magnitude_curves[i], color='red', alpha=0.7, 
                label=('' if i==0 else '_') + 'Outliers de Magnitud')
    
    # Grafica solo 3 shape outliers
    for i in range(len(shape_curves)):
        plt.plot(shape_curves[i], color='blue', alpha=0.7, 
                label=('' if i==0 else '_') + 'Outliers de Forma')
    
    plt.title('Outliers de Magnitud y Forma')
    plt.xlabel('Longitud de onda (nm)')
    plt.ylabel('Intensidad')
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()
    
    # Imprime estadísticas de outliers
    print(f"Outliers de magnitud: {len(magnitude_outliers)}")
    print(f"Outliers de forma: {len(shape_outliers)}")
    print(f"Shape outliers mostrados: {min(3, len(shape_outliers))}")

# Ejemplo de uso:

plot_kind_class_outliners(magnitude_out, shape_out, smoothed_result[:, :, wavelength_idx])

Outliers de magnitud: 3
Outliers de forma: 22
Shape outliers mostrados: 3

import numpy as np
import matplotlib.pyplot as plt

def shape_outliers(magnitude_outliers, shape_outliers, X):
    """
    Plot outliers and non-outliers with different colors and labels.
    Shows only 3 shape outliers and all other curves in grey with alpha=0.5
    
    Args:
        magnitude_outliers (list): List of indices for magnitude outliers
        shape_outliers (list): List of indices for shape outliers
        X (numpy.ndarray): Array with all curves data
    """
    # Convierte los índices de outliers a conjuntos para una búsqueda eficiente
    magnitude_outliers_set = set(magnitude_outliers)
    shape_outliers_set = set(shape_outliers)
    
    # Crear máscaras para diferentes tipos de outliers
    magnitude_outliers_mask = np.array([i in magnitude_outliers_set for i in range(len(X))])
    shape_outliers_mask = np.array([i in shape_outliers_set for i in range(len(X))])
    
    # Encuentra los magnitude outliers en X
    magnitude_curves = X[magnitude_outliers_mask, :]
    
    # Encuentra los shape outliers en X (limitado a 3)
    shape_curves = X[shape_outliers_mask, :]
    shape_curves = shape_curves[19:] if len(shape_curves) > 3 else shape_curves
    
    # Encuentra las curvas que no son outliers
    non_outliers_mask = ~(magnitude_outliers_mask | shape_outliers_mask)
    non_outliers = X[non_outliers_mask, :]
    
    # Grafica los tipos de curvas
    plt.figure(figsize=(12, 6))
    
    # Primero grafica todas las curvas no outliers en gris
    for i in range(25):
        plt.plot(non_outliers[i], color='gray', alpha=1, label='_')
    
    
    # Grafica solo 3 shape outliers
    for i in range(len(shape_curves)):
        plt.plot(shape_curves[i], color='red', alpha=1, 
                label=('' if i==0 else '_') + 'Outliers de Forma')
    
    plt.title('Outliers de Magnitud y Forma')
    plt.xlabel('Longitud de onda (nm)')
    plt.ylabel('Intensidad')
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()
    
    # Imprime estadísticas de outliers
    print(f"Outliers de magnitud: {len(magnitude_outliers)}")
    print(f"Outliers de forma: {len(shape_outliers)}")
    print(f"Shape outliers mostrados: {min(3, len(shape_outliers))}")

# Ejemplo de uso:

shape_outliers(magnitude_out, shape_out, smoothed_result[:, :, wavelength_idx])

Outliers de magnitud: 3
Outliers de forma: 22
Shape outliers mostrados: 3